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1 .  INTRODUCTION 


The  AFGL  Statistical  Analysis  Program  (ASAP)  described  by  Norquist  (1986, 
1988)  provides  an  optimal  Gandin  type  analysis  of  height  and  winds  on  the 
sigma  levels/layers  of  the  AFGL  Global  Spectral  Model  (GSM).  ASAP  follows  Dey 
and  Morone  (1965)  and  Bergman  (1978,  1979)  fairly  closely.  Detailed 
documentation  of  ASAP  are  given  in  an  earlier  report  (Hoffman,  1987). 

The  present  report  documents  changes  to  ASAP  and  experiments  designed  to 
test  the  impact  of  these  changes.  During  preparation  of  Hoffman  (1987),  we 
noted  some  possible  enhancements  and  minor  corrections  which  should  be  made  to 
ASAP.  In  the  past  year,  a  number  of  important  enhancements  have  been  made  and 
tested.  The  current  report  documents  these  activities.  Similar  changes  have 
been  made  to  the  moisture  analysis  program. 

2 .  ENHANCEMENTS 

Key  enhancements  Include  a  new,  simpler  and  more  efficient  method  of 
converting  pressure  layer  temperatures  into  sigma  level  heights,  a  corrected 
and  improved  handling  of  the  least  squares  equations  and  a  new  solver  for 
these  systems  of  equations  based  on  standard  LINPACK  routines. 

2.1  Satellite  temperature  profile  preprocessing 

Satellite  temperature  profile  data  is  generally  reported  as  layer  mean 
temperatures  or  equivalent  height  differences  between  mandatory  pressure 
levels.  These  data  must  be  anchored  to  the  forecast  height  field  (as  is  done 
in  ASAP)  or  to  a  preexisting  surface  pressure  or  1000  mb  height  analysis. 
Additionally,  these  data  must  be  interpolated  to  the  sigma  levels  of  the 
model.  Although  this  vertical  interpolation  might  be  accomplished  sta¬ 
tistically,  ASAP  interpolates  all  profile  data  analytically  in  the  vertical. 
For  satellite  temperature  profile  data  we  now  accomplish  this  transformation 
by  converting  the  mean  pressure  layer  temperature  (T^)  to  an  equivalent 
temperature  profile  defined  by  the  values  of  pressure  level  temperatures  (Ta) 
making  use  of  the  assumption  that  temperature  is  piecewise  linear  in  fn  p. 
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This  temperature  profile  is  then  integrated  hydrostatically  to  obtain  sigma 
layer  temperatures  (Tj^)  and  the  sigma  level  heights  .  In  principle,  this 
transformation  is  not  necessary,  the  analysis  scheme  could  make  use  of  the 
fact  that  is  correlated  with  T^.  In  practice,  it  is  desirable  to  make  this 
transformation,  since  the  correlations  with  individual  are  weak. 

By  consistently  making  the  assumption  that  T  is  linear  in  in  p  between 
the  p^,  we  will  introduce  no  height  errors  due  to  our  transformations  for  any 
layer  which  has  been  integrated  over.  We  note  that  in  converting  layer  to  le¬ 
vel  temperatures,  there  is  an  extra  degree  of  freedom  in  the  level  representa¬ 
tion.  We  eliminate  this  by  assuming  that  our  piecewise  linear  temperature 
profile  has  the  same  slope  in  the  two  layers  closest  to  400  mb  (p^) . 

The  subscript  notation  used  here  is  the  following;  i  is  reserved  for 
pressure,  k  for  sigma.  A  caret  over  the  first  character  of  a  subscript  indi¬ 
cates  a  level  quantity.  Subscripts  increase  with  increasing  height  (de¬ 
creasing  pressure).  Layer  i  is  bounded  by  levels  £  and  i-1.  The  first  level 
has  an  index  of  zero.  We  use  q  to  denote  in  p  throughout. 

We  assume  here  that  the  satellite  temperatures  are  actually  thicknesses 
converted  to  temperatures  hydrostatically.  That  is,  they  are  temperatures 
averaged  with  respect  to  in  p  (L.  McMillin,  pers.  comm.,  12  Feb.  88).  Since  T 
is  linear  on  q. 


T,  -  <TJ  . 


Therefore,  given  the  T^  and  one  Ta,  we  can  then  step  up 


,  “  2T,  ,  -  Ta 
i+1  i+1  i 


and  down 


Ta  ,  -  2T  -  Ta 
i-1  i  i 


to  determine  all  the  Ta.  If  p^  lies  in  layer  i,  i.e.  if  pA  ^  >  p|\  >  pA,  then 
may  calculate  T^  by  interpolating  in  q; 


we 


T|J  -  -  qj) 
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To  determine  the  first  Ta,  we  find  the  ra  such  that  p^j  minimizes  |Pjjj  -  4001 

assume  that  T  -  a  +  bq  in  layers  m  and  m  +  1.  Thus,  T^  is  the  temperature 

at  (qA  +  qA  , )/2  and  T  is  the  temperature  at  (qA  +  qA  )/2.  As  above,  we 
ra  m-1  m+1  m  m+1 

interpolate  in  q  to  obtain 


Ta  -  T  +  ((T  -  T  )/(qA  -  qA  ) ] (qA  -  qA  ) . 

m  m  m+1  m  m+1  m— 1  m  m— 1 

Given  the  Ta  and  Ta,  we  calculate  the  sigma  layer  mean  temperatures  by 

dividing  the  sigma  layer  by  the  pressure  levels.  Ue  calculate  Tj^  only  if  the 

kth  sigma  layer  is  within  our  temperature  profile,  i.e.  if  pg  >  p^  ^  . 

and  pA  >  pA.  In  general,  sigma  layer  k  will  be  covered  by  pressure  layers 
K  L 

i-m,n.  The  contribution  of  layers  m  and  n  may  be  due  to  only  part  of  the 
layer.  In  all  cases,  the  contribution  to  J  T  dq  from  any  one  layer  (or  par¬ 
tial  layer)  will  be  the  thickness  of  the  layer  times  the  average  temperature 
of  the  layer  which  is  the  average  of  the  temperature  at  the  two  end  points. 

If  B  and  T  indicate  bottom  and  top,  then 


-  ■’£>  ■ 


f“"m 


where 


1 

cr 

' - y - 

if  i  -  m 

q  -  • 

cr 

if  i  -  n 

t  qA 

otherwise . 

t  q^ 

otherwise 

f  Ta 

if  i  -  m 

H 

H 

1 

_ A _ 

1 

if  f  -  n 

H 

1 

otherwise . 

1 

L  Ta 

otherwise 

Note  that  m  might  be  equal  to  n  if  the  sigma  layer  is  within  the  pressure 
layer.  In  general  we  have: 
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k 


n 


Once  we  have  the  values,  we  obtain  the  heights  hydrostatically  by 

if  Tj^  is  not  missing,  and  Za  -  missing  if  Tj^  is  missing.  The  anchor  for  level 
k  is 

Z  if  k  -  1 

* 

Z  -  •  Z^  ^  if  Z^  ^  is  not  missing 

2a  otherwise 

k 

where  2a  Is  the  first  guess  value. 

We  have  implemented  these  suggestions  in  subroutine  ANCHOR  which  is 
called  by  a  new  version  of  subroutine  SaTLTMP  which  directly  replaces  the  old 
version.  Comparisons  of  the  old  and  new  versions  (2.0  and  3.0)  of  SATLTMP  are 
reported  below  (Section  3.1)  and  the  ANCHOR  and  SATLIMF  (version  3.0) 
procedures  are  described  in  detail  in  the  Appendix. 

2.2  Specifying  the  normal  equations 

Three  inconsistencies  in  specifying  the  normal  equations  have  been 
corrected. 

1)  The  geostrophlc  factor  acts  to  decouple  the  wind  and  height  analysis  in 
the  equatorial  belt.  A  strict  decoupling  is  now  enforced  in  this  belt. 

If  desired,  it  would  now  be  simple  to  decouple  the  wind  components  from 
each  other  to  analyze  a  divergent  wind  field.  As  of  Sept^nber,  198<i,  the 
NMC  tropical  wind  analysis  is  univariate  in  the  equatorial  belt  (Dey  et 
al . ,  1985a,  Section  D.2). 
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Previously,  Intercorrelations  poleward  of  10*  resulted  in  analyses  which 
were  not  strictly  univariate  in  Z  and  bivariate  in  (u,v).  To  force 
decoupling  some  of  the  analysis  weights  were  set  to  zero  (see  ASAP2(A0) 
in  Hoffman,  1987).  Consequently  the  estimate  of  Ag  produced  was  not 
strictly  correct,  since  the  weights  used  are  no  longer  the  least  squares 
solution.  That  is,  in  this  case  Bergman's  Eq.  2.8  and  2.11  do  not  hold 
and  his  Eq.  2.6  should  be  used  to  estimate  A£  exactly. 

For  analysis  grid  points  equatorward  of  10*  we  now  force  a  strict,  proper 
univariate  Z  analysis  and  bivariate  (u,v)  analysis  by  setting  all  uz ,  vz , 
zu  and  zv  correlations  to  zero  in  the  observation-observation  matrix  and 
the  observation-grid-point  RHS  (right  hand  sides) .  This  is  equivalent  to 
solving  two  separate  systems  of  equations  since  the  matrix  is  now  a 
permutation  of  a  block  diagonal  matrix.  The  approach  we  have  taken 
minimizes  changes  to  the  code  but  involves  some  extra  computations  (with 
zeros)  since  factorizing  a  matrix  goes  like  N**3 .  Since  the  RHS  have 
zeros  for  the  crossterms  we  are  assured  that  the  corresponding  analysis 
weights  will  be  zero.  However  we  have  the  convenience  of  retaining  the 
preexisiting  data  structure  since  all  the  data  are  kept  in  the  update. 

We  have  actually  set  the  criterion  for  equatorial  belt  as  11*  so  that  the 
Gaussian  grid  points  at  10.08*  are  included. 

2)  The  calculation  of  the  RHS  has  been  corrected.  In  the  selection  phase 
the  RHS  are  normalized  by  (1  +  t  ),  i.e.,  by  the  diagonal  of  the 
correlation-correlation  matrix.  The  normalized  RHS  are  the  analysis 
weights  we  would  obtain  if  each  observation  were  the  only  one.  However 
this  normalization  is  not  now  applied  to  the  actual  RHS  of  the  normal 
equations . 

3)  An  underground  analysis  is  performed  so  that  the  layer  Z  analysis  may  be 
Interpolated  to  level  T  increments  which  are  then  converted  to  layer  T 
increments  by  the  Flattery  algorithm. 

Although  only  a  Z  underground  analysis  is  required  we  now  perform  a 
complete  <Z,u,v)  underground  analysis.  We  want  the  lowest  two  layer  and 


5 


Che  underground  analyses  to  be  as  similar  as  possible  Co  assure 
reasonable  agreement  between  the  surface  and  first  level  temperatures. 
Therefore  our  procedures  are  exactly  the  same  for  the  underground  layer 
as  for  the  other  layers.  This  includes  data  selection  and  the  solution 
of  the  normal  equations,  including  data  rejection  as  necessary  (described 
below),  even  if  only  the  u  or  v  analysis  is  not  reasonable. 

Previously,  the  underground  analysis  was  handled  specially;  Potentially, 
the  underground  Z  analysis  could  have  been  very  different  from  the  lowest 
layer  analysis  because  of  differences  in  data  selection. 

2.3  Solving  the  normal  equations 

The  least  squares  analysis  or  normal  equations  are  solved  at  each  grid 
point  for  (Z,u,v).  It  is  desirable  that  precisely  the  same  data  are  used  for 
each  of  Z,  u  and  v  so  that  the  increments  are  geostrophically  balanced. 

Previously,  the  Choleski  algorithm  as  Implemented  by  Stobie  (1984;  Stobie 
et  al.,  1985)  was  used  to  solve  the  normal  equations.  The  version  of  this 
algorithm  we  had  been  using  does  not  properly  check  for  near  singular  or 
non-positive  definite  systems.  Consequently,  some  of  the  analysis  weights  a', 
were  not  properly  determined.  In  some  cases  these  incorrect  weights  may  be 
quite  large  in  magnitude  and  lead  to  impossible  values  of  the  estimated 
analysis  error  A£.  These  conditions  were  checked  for  and  if  necessary  an 
observation  was  deleted  and  the  normal  equations  resolved.  The  criteria  for 
deleting  an  observation  are  discussed  at  greater  length  in  Section  3.2.1. 

Since  each  of  Z,  u  and  v  were  solved  for  independently,  different  data  might 
be  used  for  each.  Worse,  improper,  but  reasonable  weights  were  used  in  some 
cases . 

The  solution  of  the  normal  equations  is  now  implemented  as  decomposition 
T 

into  R  R  form  followed  by  a  forward -backward  solution  for  each  RHS .  We  used 
the  LINPACK  routines  SPOCO  and  SPOSL  (Dongarra  et  al . ,  1979)  to  factor  and 
solve  the  system.  In  addition  the  LINPACK  routines  are  used  to  calculate 
RCOND,  an  estimate  of  the  reciprocal  of  the  condition  number  of  the  matrix. 
Calculating  RCOND  Involves  some  extra  calculations.  If  the  matrix  is  of  order 
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N  then  SPOCO  is  reported  to  take  (1  +  N/18)  times  the  time  to  simply  factor 
the  matrix. 

If  RCOND  is  too  small,  say  less  than  0.001,  then  we  eliminate  an  obser¬ 
vation  and  recalculate  the  solution.  This  is  similar  to  what  was  done  before, 
but  uses  a  more  rational  criterion.  In  choosing  the  cutoff  for  RCOND  we  have 
the  following  guidance:  If  we  are  solving  Ax  -  b  and  if  the  reciprocal  condi¬ 
tion  number  is  equal  to  10'*^  then  we  expect  the  solution  x  we  calculate  to 
have  d  fewer  significant  digits  than  the  elements  of  A.  While  we  calculate 
the  elements  of  A  to  full  machine  precision  this  calculation  is  based  on 
simple  models  of  the  statistical  functions.  It  is  difficult  to  argue  that  we 
really  know  these  correlations  to  much  more  than  2  decimal  places.  If  re¬ 
ciprocal  condition  number  is  0.001  in  a  particular  case  than  we  have  no  actual 
knowledge  of  the  proper  weights  to  use. 

Of  course  the  value  RCOND  from  SPOCO  is  an  estimate.  To  be  safe  we 
retain  loose  checks  on  a'  and  and  a  check  on  Ag.  To  avoid  updating  the 
pointers  and  creating  secondary  arrays  we  leave  a  placeholder  equation  for 
each  observation  deleted.  For  the  ith  observation  this  equation  is  a^^  -  0. 
That  is  we  set  all  correlations  with  the  deleted  observation  to  zero  and  leave 
a  one  on  the  diagonal.  Obviously,  we  will  obtain  a  zero  analysis  weight  for 
this  observation. 

As  a  possible  alternative  we  note  that  the  factorizations  produced  by  the 
LINPACK  routines  allow  one  to  solve  any  leading  subproblem.  That  is,  the 
factorizations  could  be  used  to  solve  the  system  obtained  by  dropping  the  last 
one  or  several,  rows  and  columns.  If  data  were  presorted  so  that  the  first 
datum  to  be  dropped  (if  it  proved  necessary)  were  that  last  datura  in  the 
normal  equations  and  the  second  to  be  dropped  were  the  next  to  last  datum, 
etc.,  then  only  a  single  decomposition  would  be  required.  However  this 
presorting  would  Ifkely  be  more  expensive  than  recalculating  just  a  small 
percentage  of  solutions  as  is  done  now. 


2.4  Ocher  minor  Issues 


A  number  of  other  minor  issues  have  been  dealt  with.  These  include: 

1)  The  translation  of  Quality  Indicators  (QI)  from  FGGE  Level  II  data  to 
internal  ASAP  Quality  Marks  (Q)  has  been  corrected  to  properly  account 
for  TWOS,  DROPWINDSONDE  and  COLEBA  data.  Subroutine  CQCV  performs  this 
translation.  Subroutine  CQCV  has  been  rewritten  as  described  in  the 
appendix.  Here  we  describe  our  interpretation  of  Table  IV  of  Appendix  A 
of  WMO  (1986,  pp.  5-7)  which  explains  the  QI  for  type  1  and  type  2  data. 

QI  is  composed  of  two  digits  denoted  IH  and  IV,  nominally  for  horizontal 
and  vertical  checks.  For  radiosondes,  pilot  balloons  and  aircraft  data 
this  translation  is  discussed  by  Hoffman  (1987,  pp.  81-82).  For  TWOS, 
DROPWINDSONDE  and  COLEBA  data  the  IH  digit  indicates  instrument 
uncertainty  1  being  very  low  and  9  being  very  high  with  0  indicating 
unknown,  and  the  IV  digit  indicates  a  measurement  type  with  1  being  very 
precise  and  larger  numbers  Indicating  a  larger  averaging  or  interpolation 
interval.  Again  0  indicates  unknown. 

We  adopt  the  following  translation  (see  Table  1): 

1)  For  TWOS  and  DROPWINDSONDE  we  flag  data  as  suspect  if 

IH  >  5  (uncertainty  in  V  >  8  m/s;  in  Z  >  300  m;  in  T  >  3.2 
degrees;  in  RH  >  16  %) 

or  IV  >  4  (a  vertical  interval  of  more  than  1.1  km). 

2)  For  COLEBAs  we  flag  data  as  suspect  if 

IH  >  5  (uncertainty  in  V  >  8  m/s) 
or  IV  >  1  (not  a  standard  Doppler  solution) 

3)  For  both,  an  IV  -  0  flag  is  ignored  and  the  quality  mark  is  set  on 
the  basis  of  IH.  Otherwise  we  translate  IH  and  IV  separately  and 
take  the  maximum  as  the  quality  mark.  (Note  (*)ed- items  in 
translation  tables  account  for  the  special  cases.) 

NOTE:  We  have  lumped  idsi  -  13  and  14  as  TWOS  data,  although  WMO  (1986) 
is  not  precise  about  idsi  -  13  quality  indicators. 


Table  1.  Translation  of  quality  Indicators  (QI)  into  quality  marks  (QM) 


QI 

0 

1 

2 

3 

4 

5 

6 

7 

8 

9 

idsi 

IH  QM 

*0 

0 

6 

8 

2 

5 

3 

4 

7 

99 

11,  12,  20-29 

1 

0 

0 

0 

0 

0 

6 

6 

6 

6 

13,  14,  15 

1 

0 

0 

0 

0 

0 

6 

6 

6 

6 

16 

1 

1 

1 

1 

1 

1 

1 

1 

1 

99 

all  others 

IV  QM 

1 

0 

6 

8 

2 

5 

3 

4 

7 

99 

11. 

12,  20-29 

*0 

0 

0 

0 

0 

6 

6 

6 

6 

6 

13, 

14,  15 

*0 

0 

6 

6 

6 

6 

6 

6 

6 

6 

16 

1 

1 

1 

1 

1 

1 

1 

1 

1 

99 

all 

others 

2)  U.S.A.  Special  Effort  quality  controlled  and/or  corrected  satellite 
temperature  profiles  are  now  preferred  to  raw  reports.  Previously,  the 
selection  was  the  other  way  around,  effectively  deleting  all  special 
effort  reports.  Tables  37  and  38  of  Appendix  A  of  WMO  (1986)  are  used. 
See  MASTOR^  in  the  appendix  for  details. 

3)  U.S.A.  Special  Effort  quality  controlled  cloud  drift  winds  are  now 
preferred  to  raw  reports.  Previously,  the  selection  was  the  other  way 
around.  Table  36  of  Appendix  A  of  WMO  (1986)  includes  the  codings  of  the 
USA  Special  Effort.  The  translation  used  and  other  details  are  described 
in  MAST0R6  in  the  appendix. 

4)  The  calculation  of  uu  and  w  forecast  error  correlation  used  in  the  buddy 
check  criterion  has  been  corrected.  In  FLAGS(19),  these  quantities  are 
now  calculated  without  taking  their  absolute  values. 

5)  Radiosonde  dewpoint  temperatures  are  now  converted  to  RH  in  MAST0R1(8) 
and  then  interpolated  linearly  in  In  p  in  PTOSIG(20) . 

6)  Data  rejected  in  FLAGS  by  the  gross  and  buddy  checks  are  reported  in  a 
new  file  on  unit  9  by  REJECT. 
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3 .  TESTS 


We  conducted  several  offline  tests  of  the  new  versions  of  procedures 
SATLTMP  and  ASAP2  to  verify  their  correctness,  compare  their  results  to 
results  of  the  earlier  versions  and  to  tune  certain  control  constants.  We 
then  conducted  a  one  week  data  assimilation  experiment,  using  simulated  data 
to  examine  the  impact  of  the  enhanced  ASAP. 

3.1  Testing  SATLTMP  in  GRDSTA 

In  order  to  test  the  new  version  of  SATLTMP  (version  3.0),  we  compared 
both  it  and  the  old  version  (version  2.0)  to  procedure  PTOSIG.  In  the  tests 
PTOSIG  was  presented  with  a  radiosonde  profile,  that  is  height  (Z)  and 
temperature  (T)  at  a  series  of  mandatory  (and  possibly  significant)  pressure 
levels  (p) ,  to  Interpolate  to  Z  values  on  the  sigma  level  and  T  values  in  the 
sigma  layers.  Then  each  version  of  SATLIMP  was  presented  with  a  profile  of 
pressure  layer  mean  temperatures  derived  hydrostatically  from  the  radiosonde 
heights,  to  Interpolate  to  sigma  heights  and  temperatures.  These  latter 
values  were  then  compared  to  the  corresponding  values  obtained  from  PTOSIG 
(considered  to  be  the  truth  here)  and  rrase  and  bias  statistics  accumulated. 
Before  summarizing  the  results  we  briefly  describe  the  three  procedures  and 
their  underlying  assumptions. 

Procedure  PTOSIG  was  developed  at  NMC  and  is  described  in  detail  in 
Gerlach  (1983,  pp.  44-46).  PTOSIG  assumes  temperature  (T)  is  a  linear  func¬ 
tion  of  In  p  to  interpolate  heights  from  pressure  levels  to  sigma  levels. 

Sigma  level  heights  are  then  differenced  to  give  sigma  layer  temperatures.  To 
interpolate  height  within  any  pressure  layer  PTOSIG  uses  the  heights  at  the 
bounding  interfaces  and  the  temperature  difference  across  the  layer. 

SATLTMP  version  1.0  was  based  on  the  Flattery  scheme  as  described  by 
Norquist  (1986,  pp.  55-58).  The  Flattery  scheme  finds  a  least  squares 
solution  to  two  conflicting  assumptions  about  the  temperature  profile,  namely 
that  layer  temperature  is  the  average  of  the  bounding  interface  values  and 
that  the  interface  values  may  be  obtained  by  interpolating  the  layer  values 
linearly  in  In  p. 
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SATLTMP  version  2.0  Is  an  alternative  offered  by  Norquist  (1986, 
pp.  64-65)  in  which  layer  temperatures  are  assumed  to  be  mass  weighted  layer 
averaged  temperatures  and  temperature  within  a  layer  is  assumed  to  vary 
linearly  with  In  p.  However  as  noted  above,  satellite  profile  layer 
temperatures  are  really  height  differences  expressed  as  (virtual)  temperatures 
hydrostatically,  i.e.  averaged  with  respect  to  In  p,  not  p.  Version  2.0  still 
uses  the  Flattery  scheme  to  get  one  sigma  level  temperature  near  400  mb. 

SATLTMP  version  3.0,  described  above  in  Section  2.1,  assumes  that 
temperature  varies  linearly  with  In  p  within  each  pressure  layer.  It  is  also 
assumed  that  this  relationship  is  differentiable  across  the  pressure  level 
closest  to  400  mb. 

In  the  following  test  cases  both  PTOSIG  and  SATLTMP  were  presented  data 
defined  at  or  between  the  mandatory  pressure  levels  1000,  850,  700,  500,  400, 
300,  200,  100,  70,  50.  These  tests  were  conducted  by  making  minor  changes  to 
the  grid  to  station  verification  program  of  Norquist  (pers.  comm.,  1987). 
Results  of  the  test  of  SATLTMP  versions  2.0  and  3.0  (SATONE  and  SATTWO)  are 
compared  here.  The  sample  runs  from  79  February  8  1200  GMT  through  February 
15  0000  GMT  excluding  February  13  1200  GMT. 

The  numerical  results  are  displayed  in  Tables  2  and  3  for  height  at  sigma 
levels  and  temperature  at  sigma  layers  respectively.  The  results  are  also 
displayed  graphically  in  Figs.  1  and  2.  Clearly,  except  for  the  eleventh 
level/layer,  both  schemes  are  adequate  and  roughly  comparable.  However 
version  3.0  has  smaller  biases  and  rmse  and  is  definitely  better  for  the 
highest  level/layer.  Version  3.0  Is  also  simpler  in  its  assumptions  and 
coding,  and  more  efficient  as  judged  by  the  computer  timings  of  the  test 
cases.  Timings  for  a  single  synoptic  time  give  2.79  and  2.43  cp  seconds  of 
execute  time  for  versions  2.0  and  3.0  respectively  for  a  net  savings  of  more 
than  10%.  However  version  3.0  does  not  process  as  many  levels  as  version  2.0. 

3.2  Testing  new  CHLSKY  and  ASAP2 

3.2.1  Analysis  of  recalculation  criteria 

To  see  the  relationship  between  old  and  new  criteria  for  rejecting  data 
and  recalculating  the  weights,  we  have  extracted  all  data  needed  to  rerun  the 
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Table  2.  Height  differences  statistics  (m)  comparing  SATLTMP  Versions  2.0 
and  3.0  versus  PTOSIG. 


LEVEL  SIGMA 

Version  2.0 

Version  3.0 

RMSE  BIAS 

SAMPLE 

RMSE 

BIAS  SAMPLE 

1 

.925 

2.9253 

-.0917 

2270 

2.8628 

.2145 

2245 

2 

.800 

9.4633 

2.8972 

4478 

5.4872 

-  .2177 

4411 

3 

.650 

9.2876 

2.1642 

7152 

5.7740 

.6397 

7015 

4 

.500 

8.7859 

.9099 

7484 

4.2194 

-  .0470 

7351 

5 

.375 

12.7985 

3.0962 

7827 

10.2658 

1.6389 

7815 

6 

.300 

14.9905 

1.1486 

7719 

13.5419 

1.9475 

7677 

7 

.250 

18.9387 

- .8738 

7753 

17.2810 

2.4887 

7703 

8 

.200 

23.6651 

-6.7386 

7185 

20.2309 

-1.6692 

7139 

9 

.150 

35.8403 

-14.3255 

7054 

30.8529 

-7.0962 

7002 

10 

.100 

46.9694 

-5.9414 

6161 

16.7813 

2.3607 

6110 

11 

.050 

123.2643 

-43.2206 

1574 

8.8502 

-  .8772 

1445 

12 

.010 

- 

- 

0 

- 

- 

0 

Table  3.  Temperature  difference  statistics  (K)  comparing  SATLTMP  Versions 


2 . 0  and 

3.0  versus  PTOSIG. 

LEVEL  SIGMA 

Version  2.0 

Version  3.0 

RMSE  BIAS 

SAMPLE 

RMSE 

BIAS  SAMPLE 

1 

.962 

1.2810 

-  .0402 

2270 

1.2537 

.0939 

2245 

2 

.861 

2.0373 

.2547 

4478 

1.4436 

- .0742 

4411 

3 

.724 

1.4712 

-  .0742 

7152 

1.1248 

.1265 

7015 

4 

.573 

.6891 

-  .1486 

7484 

.5886 

- .0897 

7351 

5 

.436 

.6426 

.2511 

7827 

.7287 

.1679 

7815 

6 

.337 

1.0466 

-  .2401 

7719 

1.1015 

.0632 

7677 

7 

.274 

1.5945 

-  .2624 

7753 

1.6635 

.0996 

7703 

8 

.224 

2.6612 

-  .9474 

7185 

2.3219 

-  .6406 

7139 

9 

.174 

2.7856 

-  .8259 

7054 

2.6414 

-  .6442 

7002 

10 

.124 

4.6059 

.6537 

6161 

2.8510 

.7667 

6110 

11 

.730 

6.1634 

-1.6887 

1574 

.6359 

- .0787 

1445 

12 

.020 

- 

0 

- 

- 

0 

Height  differences  statistics  (m)  comparing  SATLTMP  Versions  2.0  and  3.0.  Truth  is  take 
radiosondes  interpolated  by  PTOSIG.  The  sample  consists  of  one  week  from  SOP-1. 


least  squares  problem  for  a  subset  of  the  analyses  for  one  synoptic  time  (06 
GMT  17  June  1979)  from  our  STATSAT  experiment  (Louis  et  al . ,  1987).  The  sub¬ 
set  in  question  includes  all  systems  of  equations  for  a  single  level  (K  -  5) 
which  did  not  yield  satisfactory  solutions  according  to  the  old  criteria.  Z, 
u  and  V  problems  were  included.  At  level  5,  out  of  a  total  of  61x62  grid 
points  there  were  98  problem  grid  points.  There  were  a  total  of  281  least 
square  systems  in  our  subset.  For  each  problem  grid  point  two  or  all  three  of 
the  systems  for  2,  u  and  v  may  be  present  and  there  may  be  multiple  systems 
for  each  if  it  was  necessary  to  remove  more  than  one  datum  before  a  satis¬ 
factory  solution  was  obtained. 

For  these  cases,  we  have  examined  certain  key  quantities  and  their 
inter-relationships.  These  quantities  are; 

Ag  the  normalized  analysis  error 

PIVOT  the  least  pivot  used  by  the  CHLSKY  subroutine 

AP  the  largest  absolute  value  of  the  normalized  weights 

RCOND  the  reciprocal  of  the  condition  number  calculated  by  LINPACK 
routine  SPOCO.  We  set  RCOND  to  0  if  the  matrix  is  singular. 

(Note;  most  of  the  plots  identify  RCOND  as  condition  number.) 

The  original  ASAP2  determines  a  solution  is  unsatisfactory  if  1)  CHLSKY  finds 
a  pivot  whose  absolute  value  is  less  than  10*^;  or  2)  AP  is  larger  than  1.1; 
or  3)  A£  Is  not  within  [0,1].  In  fact  condition  (1)  never  occurred  in  our 
sample,  although  there  were  many  cases  of  negative  PIVOT.  Choleski  pivots 
should  never  be  zero  or  negative.  This  indicates  a  nonpositive  definite 
matrix.  Criterion  (2)  is  arbitrary;  some  solutions  with  AP  of  1.5  or  2.0 
might  be  valid,  while  some  solutions  with  small  values  of  AP  are  actually 
poorly  conditioned.  Condition  (3)  only  occurred  when  (2)  had  also  occurred  in 
the  cases  studied. 

We  note  from  the  plots  in  Fig.  3  that  large  values  of  AP  are  associated 
with  pivots  of  small  absolute  magnitude.  Large  values  of  AE  are  associated 
with  small  negative  values  of  PIVOT  and  with  large  values  of  AP.  We  found 
that  most  of  the  cases  of  negative  and  small  Ag  (less  than  say  0.2)  are 
associated  with  10.08*N,  53.11'E.  These  points  are  circled  in  some  of  the 
plots.  We  were  unable  to  determine  the  difficulty  with  these  systems. 


:atterplots  desc 
riginal  version 
>ry.  See  text  f' 


We  see  from  the  plots  in  Fig.  U,  that  condition  number  is  closely  related 
to  PIVOT  and  that  large  values  of  AP  and  A£  generally  have  small  RCOND.  Using 
a  cutoff  of  0.01  for  RCOND  would  trap  all  but  a  few  of  the  cases  trapped  by 
the  old  scheme  in  this  sample. 

We  find  that  each  time  RECALC,  the  procedure  to  remove  one  of  the  most 
highly  correlated  data  in  the  system,  is  invoked  the  condition  number  of  the 
matrix  improves,  the  AP  maximum  value  decreases  and  the  A^  value  increases. 

The  magnitude  of  these  changes  seems  correlated.  That  is,  if  dropping  a  point 
significantly  improves  the  condition  number,  it  also  makes  a  large  difference 
to  AP  and  AE. 

Since  the  relationship  between  RCOND  and  AP  is  not  well  defined  it  seems 
preferable  to  eliminate  data  on  the  basis  of  RCOND  which  has  a  well  defined 
mathematical  meaning  even  if  the  cutoff  value  is  chosen  empirically.  Fig.  5 
shows  that  the  estimate  of  reciprocal  condition  number  from  SPOCO  compares 
well  with  the  estimate  from  SCHDC,  another  LINPACK  routine  after  we  account 
for  the  fact  that  the  estimates  from  SCHDC  are  generally  overestimates,  often 
by  a  factor  of  10  (Dongarra  et  al.,  1979,  p.  8.14). 

3.2.2  Tuning  of  new  version  of  CHLSKY 

The  new  version  of  CHLSKY  has  two  adjustable  parameters  in  the  specifi¬ 
cation  of  the  data  rejection  criteria.  These  are  RCRIT  and  XCRIT.  When 
either  RCOND/RCRIT  is  less  than  one  (i.e.  the  system  is  poorly  conditioned)  or 
XMAX,  the  largest  normalized  data  weight,  is  greater  than  XCRIT,  the  most 
highly  correlated  datura  is  eliminated  from  the  system  and  the  reduced  system 
is  resolved.  We  conducted  several  test  runs  for  the  case  of  06  GMT  17  June 
1979  using  the  new  versions  of  ASAP2  and  CHLSKY  as  described  in  the  appen¬ 
dix.  The  statistics  reported  below  are  based  on  all  systems  of  equations  for 
this  case . 

In  the  first  run  RCRIT  -  0.001  and  XCRIT  -  4.0.  In  this  run,  except  for 
extreme  values  of  a' ,  all  data  rejection  is  based  on  the  condition  number  of 
the  system  of  equations.  While  there  were  many  cases  of  Ag  -  1 ,  i.e.  cases  of 
no  reduction  in  analysis  error,  there  were  only  4  cases  of  negative  A£.  These 
impossible  values  of  Aj.  were  associated  with  marginal  condition  numbers  (RCOND 
<  3/1000)  and  large  AP  values  (XMAX  >  3,0;  three  of  these  four  cases  had  XMAX 
>  4.0).  There  were  a  total  of  16  cases  of  XMAX  >  4.0. 
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Scatterplots  describing  relationships  among  parameters  used  by  the 
original  version  of  CHLSKY  to  check  if  the  solution  is  satisfactory 
and  reciprocal  condition  number  calculated  by  LINPACK  routine 
SPOCO.  See  text  for  details. 


SPOCO 


ScaCterplot  describing  relationship  between  reciprocal  condition 
numbers  calculated  by  LINPACK  routines  SPOCO  and  SCHDC.  See  text 
for  details. 


Out  of  a  fflaximuffl  possible  A9166  grid  points,  roughly  37000  had  some 
data.  In  total  14786  pieces  of  data  at  roughly  8100  grid  points  were 
eliminated  in  this  run.  Of  these  more  than  half,  8088  pieces  of  data  at  4440 
grid  points,  corresponded  to  non-positive  definite  systems,  i.e.  RCOND  -  0.0. 

We  repeated  the  experiment  outlined  above  with  XCRIT  -  1  and  RCRIT  - 
.00001.  These  values  cause  data  rejection  if  the  system  of  equations  is  non¬ 
positive  definite  but  not  if  it  is  ill-conditioned  unless  a  value  of  a'  >  1  is 
determined.  By  examining  the  distribution  of  XMAX  for  poor  and  well  condi¬ 
tioned  systems  we  see  that  many  poorly  conditioned  systems  have  fairly  reason¬ 
able  values  of  XMAX.  This  is  seen  in  Table  4,  which  is  a  table  of  histograms 
of  the  number  of  occurrences  of  different  values  of  RCOND  for  different  suben¬ 
sembles  defined  by  XMAX  >0,  1,  2,  3  or  4.  The  singular  cases  are  not 
included  here.  The  important  result  is  that  there  are  many  poorly  conditioned 
systems  which  have  reasonable  values  of  XMAX. 


Table  4.  Distribution  in  terms  of  reciprocal  condition  number,  conditioned  by 
size  of  maximum  normalized  weight. 


1000*RCONP _ 

Occurrences 

for  XMAX  > 

> 

<- 

0(*) 

1 

2 

3 

4 

0.0 

0.001 

74 

0 

0 

0 

0 

0.001 

0.01 

593 

1 

1 

1 

1 

0.01 

0.1 

598 

265 

218 

173 

155 

0.1 

1 

9270 

1870 

516 

336 

234 

1 

10 

24274 

2370 

234 

58 

14 

10 

100 

6402 

408 

0 

0 

0 

100 

1000 

1370 

1 

0 

0 

0 

0 

1000 

42581 

4915 

969 

568 

404 

(*)  Estimate  based  on  first  2000  systems  of  equations. 
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From  this  table  we  conclude  that  values  of  RCRIT  of  .001  and  XCRIT  of  2.0 
are  reasonable.  With  these  values  almost  all  data  rejections  will  result  from 
detecting  poorly  conditioned  systems;  only  a  few  hundred  additional  data  will 
be  rejected  by  the  other  criteria.  It  does  not  seem  worth  while  to  put  fur¬ 
ther  effort  into  optimizing  the  rejection  criteria  for  the  following  reasons: 

1)  The  data  being  rejected  are  highly  correlated  with  other  data  which  is 
kept.  It  would  be  more  worthwhile  working  on  a  methodology  which 
combines  the  two  data  rather  than  rejecting  one  (e.g.  Lorenc ,  1981). 

2)  The  majority  of  the  data  rejected  are  due  to  singular  systems,  not 
ill-conditioned  ones.  These  data  must  be  rejected  in  the  current 
formulation  of  the  01. 

3)  A  reformulation  of  the  horizontal  correlation  function  in  terras  of  Bessel 
functions  may  result  in  better  conditioned  systems,  reducing  the 
importance  of  the  rejection  criteria. 

3.2.3  Estimated  analysis  errors  and  corrections:  impact  and  sensitivity 

In  order  to  assess  the  impact  and  sensitivity  of  the  new  choleski  scheme 
on  ASAP,  we  compared  the  estimated  analysis  errors  and  corrections  from  the 
old  scheme  (CONTROL) ,  the  new  scheme  (NEW)  and  one  of  the  test  cases  (TEST) 
which  was  described  above.  NEW  incorporates  all  the  changes  discussed  in 
Section  2  (except  for  change  5  in  Section  2.4)  and  makes  use  of  RCRIT  -  0.001 
and  XCRIT  -  2.  TEST  incorporates  only  the  changes  discussed  in  Sections  2.2 
and  2.3  and  makes  use  of  RCRIT  -  0.001  and  XCRIT  -  4.  We  identify  NEW  - 
CONTROL  as  IMPACT,  since  it  is  a  measure  of  the  impact  of  all  changes  and  NEW 
-  TEST  as  SENSITIVITY  since  it  is  a  measure  of  the  sensitivity  to  small 
differences  in  data  selection  and  in  the  parameters  of  the  new  choleski 
solver.  Summary  statistics  are  displayed  in  Tables  5,  6,  7  and  8  for  height 
and  u-wind  component  only.  Results  for  v-wind  component  are  similar  to  the 
u-wind  component  results. 

Considering  the  mean  EAE  impact  we  see  the  new  scheme  results  in  esti¬ 
mates  of  the  analysis  error  which  are  smaller  by  order  one  decameter  in  height 
and  one  m/s  in  wind  components.  The  size  of  the  variation  of  IMPACT  is  the 
same  order  of  magnitude.  By  comparison  EAE  sensitivity  is  relatively  small. 
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Table  5.  Mean  heights  (m) 


Level 

EAE 

Correction 

NEW 

IMPACT  SENSITIVITY 

NEW 

IMPACT  SENSITIVITY 

1 

21.825 

-1.045 

.991 

.025 

-.192 

-  .405 

2 

24.689 

-2.294 

.107 

-.364 

-  .724 

-1.035 

3 

32.730 

-3.667 

-  .085 

-1.303 

-1.597 

1.673 

k 

37.558 

-3.932 

-  .042 

-2.034 

-1.867 

-1.879 

5 

48.215 

-4.696 

-  .058 

-2.931 

-2.584 

-2.528 

6 

61.583 

-6.258 

-.057 

-4.083 

-3.470 

-3.434 

7 

69.465 

-7.641 

-.069 

-4.417 

-3.552 

-3.615 

8 

72.536 

-8.431 

-.088 

-4.543 

-3.529 

-3.592 

9 

72.002 

-8.631 

-.122 

-5.305 

-3.572 

-3.414 

10 

73.266 

-9.135 

-.147 

-5.638 

-3.480 

-3.313 

11 

83.717 

-6.313 

-  .074 

-4.861 

-2.312 

-1.896 

12 

98.942 

-1.248 

.020 

-1.323 

.316 

.769 

Table  6. 

Standard 

deviation  of  height  (m) 

Level 

EAE 

Correction 

NEW 

IMPACT  SENSITIVITY  NEW 

IMPACT  SENSITIVITY 

1 

13.887 

1.702 

1.996 

5.402 

4.308 

4.038 

2 

16.086 

2.616 

.789 

6.810 

4.322 

3.556 

3 

19.709 

4.340 

.810 

10.240 

5.162 

4.349 

4 

23.178 

4.516 

.892 

11.931 

4.708 

3.866 

5 

28.119 

5.244 

1.098 

15.307 

5.398 

4.527 

6 

34.450 

7.089 

1.407 

19.633 

7.203 

5.910 

7 

36.519 

8.479 

1.635 

22.919 

8.308 

6.415 

8 

35.119 

9.026 

1.743 

24.735 

8.832 

6.688 

9 

31.917 

9.322 

1.832 

27.512 

9.902 

7.831 

10 

26.856 

9.403 

2.102 

31.446 

12.126 

9.237 

11 

18.541 

6.405 

1.063 

28.169 

11.151 

10.221 

12 

11.184 

2.246 

.634 

17.735 

10.509 

11.956 
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Table  7.  Mean  u-wind  component  (m/s) 


Level 

EAE 

Correction 

NEW 

IMPACT  SENSITIVITY 

NEV 

IMPACT  SENSITIVITY 

1 

4.548 

-  .145 

.058 

-.008 

.022 

.021 

2 

5.197 

-.248 

-.005 

-.027 

.006 

.019 

3 

6.982 

-.372 

-.008 

-.040 

.016 

.022 

4 

7.938 

-  .405 

-.002 

-.050 

-  .004 

.004 

5 

10.230 

-  .478 

-.002 

.009 

.011 

.019 

6 

13.072 

-  .627 

-  .002 

-.020 

.001 

.017 

7 

14.611 

-  .806 

-  .007 

-.008 

-  .001 

.028 

8 

15.307 

-.884 

-  .008 

.027 

.007 

.023 

9 

15.532 

-  .835 

-  .014 

.057 

.040 

.044 

10 

13.979 

-  .737 

-.015 

.114 

.085 

.072 

11 

15.042 

-  .559 

-  .023 

.062 

.084 

.083 

12 

13.612 

-  .189 

-.012 

.127 

-  .003 

-  .034 

Table  8. 

Standard 

deviation  of  u-wind 

component 

(m/s) 

Level 

EAE 

Correction 

NEW 

IMPACT  SENSITIVITY 

NEW 

IMPACT  SENSITIVITY 

1 

2.288 

.271 

.184 

1.170 

.690 

.569 

2 

2.695 

.343 

.154 

1.312 

.712 

.549 

3 

3.769 

.476 

.166 

1.664 

.724 

.572 

4 

4.167 

.527 

.171 

1.783 

.  644 

.497 

5 

5.336 

.612 

.206 

2.159 

.726 

.598 

6 

7.008 

.762 

.261 

2.785 

.874 

.652 

7 

7.750 

.949 

.299 

3.279 

.956 

.684 

8 

7.843 

1.007 

.296 

3.390 

1.130 

.876 

9 

7.596 

.981 

.283 

3.249 

1.102 

.825 

10 

6.373 

.905 

.251 

2.942 

1.040 

.781 

11 

6.462 

.816 

.330 

2.549 

1.052 

.975 

12 

5.281 

.579 

.265 

1.282 

.746 

.811 
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Note  that  the  changes  to  SATLTMP  effect  the  data  values  but  not  the 
observing  error  statistics.  Therefore  CONTROL  and  TEST  corrections  should  be 
quite  similar.  This  is  especially  evident  In  the  mean  height  correction 
results  where  IMPACT  and  SENSITIVITY  are  nearly  the  same.  The  height  correc¬ 
tions  themselves  are  roughly  half  the  size  of  typical  radiosonde  observing 
errors.  The  height  IMPACT  and  SENSITIVITY  are  of  order  one  decameter.  Mean 
wind  corrections  are  quite  small;  the  size  of  the  wind  IMPACT  and  SENSITIVITY 
is  of  order  one  m/s . 

We  also  examined  horizontal  maps  of  EAE  for  NEW,  IMPACT  and  SENSITIVITY 
at  .574  sigma.  IMPACT  is  small,  but  persistent  in  the  data  rich  areas  and 
substantial  and  persistent  in  the  oceanic  data  voids.  SENSITIVITY  is 
generally  quite  small  and  isolated.  SENSITIVITY  may  be  quite  large  in  small 
areas  in  equatorial  areas  such  as  the  Amazon  Basin. 

Actual  impacts  after  spectral  filtering  and  postprocessing  are  generally 
small  (see  Fig.  6).  There  are  occasional  noticeable  small  scale  features. 
However  changes  will  grow  during  the  first  few  days  of  the  assimilation  cycle 
as  the  EAE  stabilize.  That  is  the  new  EAE  will  be  used  to  forecast  the  EPE 
during  the  subsequent  analysis.  Note  that  the  EAE  is  spectrally  filtered 
during  this  process. 

3.3  Impact  test 

We  conducted  a  one  week  simulated  data  assimilation  experiment  with  the 
old  and  new  ASAP  to  study  the  impact  of  the  enhancements  described  in  Sec¬ 
tion  2.  The  control  experiment,  denoted  STATSAT,  was  previously  run  as  one  of 
the  baseline  experiments  in  our  studies  of  the  impact  of  SSM/T-l  and  SSM/T-2 
and  WINDSAT  lidar  wind  data.  The  experiment  using  the  new  ASAP  is  denoted 
OITEST.  (Change  number  5  of  section  2.4  was  not  implemented  in  OITEST,  and 
the  original  version  of  the  moisture  01  was  used  in  the  assimilation.)  The 
only  differences  between  the  two  experiments  are  the  enhancements  to  the  ASAP 
described  here.  The  nature  run  for  these  experiments  is  the  ECMWF  N48  grid 
point  model  simulation  of  10  -  30  November  1979.  The  experiments  begin  at 
0000  GMT  on  18  November.  To  obtain  realistic  initial  conditions  we  first  ran 
a  96  h  forecast  from  the  nature  run  at  0000  GMT  11  November  and  then  a  72  h 
spinup  assimilation  cycle.  During  the  assimilations  described  here  the  full 
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FGGE  Level  II  data  set  created  by  Dey  et  al .  (1985b)  exclusive  of  the  Ildar 
winds,  was  used  in  exactly  the  same  fashion  as  described  by  Louis  et  al . 

(1987). 

The  experimental  setup  used  is  ideal  for  detecting  impacts  due  to  changes 
in  the  assimilation  system.  First  of  all,  for  verification  purposes  we  know 
nature  and  can  perform  grid  point  to  grid  point  comparisons.  Second  the  data 
errors  are  ideal  and  known  exactly.  This  eliminates  the  effect  present  in 
real  data  experiments  of  not  knowing  the  errors  statistics  exactly.  On  the 
other  hand  these  idealized  errors  are  too  easy  to  assimilate  and  results  from 
such  experiments  are  overly  optimistic.  As  described  by  Dey  et  al.  (1985b, 
Table  1) ,  the  simulated  observational  errors  are  all  pure  random  numbers 
except  that  the  satellite  temperature  errors  include  a  small  level  by  level 
bias  depending  on  the  total  cloud  amount  as  estimated  from  the  relative 
humidity  profiles.  Also  the  nature  run  itself  is  smoother  than  the  real 
atmosphere.  Consequently,  the  analysis  is  able  to  average  out  the 
observational  errors  very  effectively. 

In  general  the  results  of  the  STATSAT  and  OITEST  are  very  similar.  It  is 
difficult  to  determine  which  is  better  by  examining  synoptic  fields.  For 
example  in  the  synoptic  1000  mb  and  500  mb  height  fields  at  the  end  of  the 
experiments  at  0000  GMT  25  November  it  is  evident  that  STATSAT  and  OITEST  both 
agree  very  well  with  NATURE  in  large  scale  features.  Smaller  scale  features 
tend  to  differ  somewhat  between  STATSAT  and  OITEST,  while  NATURE  has  no 
smaller  scales  at  all.  As  an  example  we  show  the  500  mb  height  fields  and 
differences  in  the  Southern  Hemisphere  in  Fig.  7.  The  wind  fields  at  200  mb 
all  look  alike.  Wind  field  errors  are  small  and  distributed  fairly  uniformly 
over  the  globe . 

To  quantify  our  results  we  calculated  global  mean  and  rms  error  statis¬ 
tics  for  STATSAT  and  OITEST.  While  OITEST  is  generally  better  than  STATSAT, 
differences  between  the  errors  are  small.  In  general  the  OITEST  wind  analyses 
show  a  small  consistent  improvement  at  all  levels.  The  OITEST  height  analyses 
are  slightly  better  at  1000  mb  but  slightly  worse  or  no  different  from  850  mb 
through  150  mb.  At  100  mb  the  OITEST  height  analysis  is  noticeably  worse.  At 
the  uppermost  two  levels,  70  and  50  mb,  the  OITEST  height  and  wind  analyses 
are  noticeably  improved.  Selected  time  evolutions  of  the  global  rms  errors 
are  shown  in  Fig.  8  and  the  time  global  rms  errors  are  displayed  in  Fig.  9. 


Southern  Hemisphere  500  mb 
hel^t  fields  for  1979  November 
17  0000  GMT:  (a)  NATURE,  (b) 
STATSAT  analysis,  (c)  OITEST 
analysis,  (d)  STATSAT  error  and 
(e)  OITEST  error. 


Fig.  8.  Time  evolution  of  global  analysis  error  during  the  STATSAT  (dashed) 
and  OITEST  (solid)  assimilation  experiments  for  rms  vector  wind 
error  (m/s)  at  (a)  1000  mb,  (b)  200  mb,  and  for  rms  height  error  (m) 
at  (c)  1000  mb,  (d)  500  mb,  (e)  100  mb  and  (f)  50  mb. 
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While  it  is  not  possible  to  determine  which  enhancements  were  most 
significant  in  this  impact  from  the  experiments  conducted,  we  hypothesize  that 
the  impact  at  the  uppermost  levels  is  due  to  the  enhanced  treatment  of 
satellite  temperature  profile  data. 

4.  CONCLUSIONS 

The  changes  made  to  ASAP  have  enhanced  its  consistency  and  operation. 
Separate  studies  of  the  satellite  profile  data  preprocessing  and  the  specifi¬ 
cation  and  solution  of  the  normal  equations  demonstrated  that  the  changes  made 
are  desirable.  The  overall  positive  impact  of  the  enhancements  is  small  as 
measured  by  a  one  week  data  assimilation  experiment  based  on  simulated  data. 

Specific  conclusions  which  may  be  drawn  from  the  experiments  described  in 
this  report  are: 

*  The  new  satellite  temperature  profile  preprocessor  SATLTMP  3.0  is 
simpler,  more  consistent  with  procedures  of  the  data  producer,  more  efficient 
in  terms  of  computer  time  and  more  accurate.  SATLTMP  3.0  is  based  on  the 
fewest  assumptions  and  is  designed  so  that  errors  made  in  individual  layers 
tend  to  cancel . 

*  The  new  specification  and  solution  of  the  normal  equations  is  more 
consistent  and  results  in  smaller  estimated  analysis  errors. 

*  The  old  method  of  solution,  while  correct  for  positive  definite 
systems,  fails  to  properly  trap  non-positive  definite  systems.  This  yields 
incorrect  solutions,  some  of  which  are  not  trapped  by  later  checks.  It  may  be 
that  these  solutions  introduced  some  noise  into  the  previous  analyses. 

*  The  one  week  impact  experiment  showed  a  significant  but  small 
improvement  due  to  the  enhancements  made.  It  is  not  possible  to  determine 
which  enhancements  were  most  significant  in  this  impact  from  the  experiments 
conducted. 
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APPENDIX.  PROGRAM  DOCUMENTATION. 


This  appendix  documents  the  new  subroutines  and  modifications  to  the  existing 
subroutines  in  the  sane  format  as  Hoffman  (1987).  Each  page  has  two  page 
numbers,  one  is  the  page  number  within  this  appendix  (e.g.  A- 3)  and  one  is  the 
insert  page  number  within  Hoffman  (1987)  (e.g.  Rev.  1.10,  21.3).  Pages  21-88 
of  D  should  be  replaced  by  the  following  pages.  Changes  and  additions  are 
marked  in  the  left  hand  margin. 
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This  page  intentionally  left  blank. 


3.  Program  Units 


The  main  program  ASAPl  and  the  various  subprograms  it  involves  are 
described  here  in  alphabetical  order.  Within  each  description,  there  are 
various  numbered  subsections.  The  start  of  each  of  these  blocks  is  identified 
by  a  comment  card  in  the  FORTRAN  code. 
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A-3 


ANCHOR 


Name : 

Purpose : 


Author : 

Documentation; 

Referenced  by: 

References : 

Commons  used: 


Arguments : 


ANCHOR  (subroutine) 

Transforms  SATEM  thickness  temperatures  into  heights.  Thick 
nesses  are  anchored  to  Z.*.  or  first  guess  height.  To  convert 
pressure  layer  thickness  to  sigma  layer  thickness,  a  piece - 
wise  linear  in  in(p)  -  q  temperature  profile  is  constructed. 

Ross  Hoffman,  AER,  1988 

Ross  Hoffman,  AER,  1988 


SATLTMP 

None 

None 

i  -  0, 

L 

in  pA,  log  pressure  at  pressure  layer 

interface  (mb) 

Ti 

i  -  1. 

L 

pressure  layer  temperatures  (K) 

L 

number  of  pressure  layers 

k  -  0, 

K 

in  p^,  log  pressure  at  sigma  layer 

interface  (mb) 

z* 

surface  height  (m) 

Za 

k 

k  -  1. 

K 

interface  heights  (m) 

k  -  1. 

K 

first  guess  interface  heights  (m) 

K 

number  of  sigma  layers 

VMISS 

missing  value 

I  FAIL 

set  to  1  if  ANCHOR  fails,  0  otherwise 

i  -  0. 

L 

workspace  for  pressure  level  temperatures 

Tk 

k  -  1, 

K 

workspace  for  sigma  layer  temperatures 

TC 

k  -  0, 

K 

workspace  for  sigma  level  temperatures 

A- 
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ANCHOR 


I/O  Units:  None 

Description: 

(0)  Parameters 

g  gravity  9.81  m/s 

R  gas  constant  287.05  (m/s)^ 

in  anchor  pressure  in  400  mb 

(1)  If  any  -  VMISS  go  to  (FAIL)ure 

(2)  Determine  m  such  that  |qA  -  q  |  is  minimized. 

m  a 

If  m  is  not  found  or  m  <  1  or  m  >  L-1  go  to  (FAIL)ure 


(3)  Determine  Ta  -  T  +  ( (T  —  T  )/(qA  ,  ~  qA  )J(qA  ~  qA  ) 
m  m  m+1  m  m+1  nr-l  m  m-1 


(4)  For  i  -  m+l,  L  determine  Ta  -  2T.  -  Ta  , 

i  i  i-1 


For  i  -  ra,  1  determine  Ta  -  2T  —  Ta 

i-1  i  i 

(5)  Determine  T/\.  First  find  layer  including  Tjv,  then  interpolate  as 
follows.  Set  i-1.  Then  for  k  —  0,  K.  If  qA  >  qA  >  qA,  then 
(else  Tj\  -  missing) 

Until  ^1  increment  i 

T^  -  Ta  +  ((TA_^  -  TA)/(q^_^  -  qj)l(q^  -  qj) 

(6)  Begin  determination  of  T  .  Find  pressure  layers  i  -  m,  n  which  cover  the 
sigma  layer,  as  follows.  Set  n  -  1,  then  for  k  -  1,  K. 


If  qA  >  qA  and  qA  >  qA,  then  (else  T  -  missing) 

Ic  ^  iC  IC 


Set  ro  -  n;  until  (qA  i]  qA  ,  >  qA  increment  ra 
nt-l  k-1  m 


Set  n  -  m;  until  (qA  >]  qA  >  qA  increment  n 
n-1  k  n 
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(7)  Accumulate  layer  contribution.  Set  -  0;  then  for  i  -  m 


If  J  -  then  <1^  -  -  Ta_j 


else  q  -  qA  ,  T  -  Ta 

^i-1  B  i-1 


If  i  -  n,  then  q  -  qA  ,  T  -  Ta 

k  T  k 


else  q  -  qA  ,  T  -  Ta 
T  i  T  t 


T  -T  +(q  -q)(T  +T) 

k  k  ^B  T  B  T 


(8)  Normalize 

T^  -  \/<2(qA_^  -  q^)):  end  k  (6). 

(9)  Calculate  Z^.  For  k  -  1,  K.  If  is  not  missing  (else 

Determine  anchor 

Z  if  k-1 
* 


Z  - 


If  Za_^  VMISS 


otherwise. 


Integrate  hydrostatically 


10)  Return  and  error  handling 


.  n 


-  missing) 
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Name : 

Purpose : 


Author: 

Documentation: 

Referenced  by: 

References : 

Commons  used: 

Arguments : 

I/O  units: 
Description: 

(1)  Initialize 

Rev  1.10 


ASAPl  (main  program) 

Performs  the  statistical  analysis  procedure  following  Bergman 
(1979)  and  Dey  and  Morone  (1985)  as  described  in  Norquist 
(1982,  1983,  1984,  1986a.  1986b)  and  Halberstam  et  al . 

(1984).  The  moisture  analysis  is  now  separate  but  moisture' 
residuals  are  calculated  here  for  the  GWE  Level  II  data  and 
stored  with  their  associated  identification  and  location 
parameters  on  units  8  and  3.  ASAPl  calls  several  routines  to 
read  in  the  data,  interpolate  to  sigma  coordinates,  and 
perform  quality  control.  It  then  loops  over  all  the  analysis 
grid  points,  selecting  observations  and  calling  ASAP2  to 
actually  perform  the  analysis.  Note  then  an  observation  may 
be  an  entire  RAOB  report  or  a  single  level  AIREP.  The  RAOB 
contains  many  pieces  of  data .  the  AIREP  only  one. 

D.  Norquist.  SASC,  1980  -  1986 

R.  Hoffman,  AER,  1986 

None 

ASAP2,  CNTOBS,  FLAGS,  MASTORl,  MAST0R2,  MAST0R4,  MAST0R6 , 
POINTS,  SETFG. 

AS140,  ASIAS2,  CCONST,  DCONST,  RESID,  SIGK 

None 

2.  3,  7,  8,  10,  INPUT,  OUTPUT 


constants.  Print  version  to  OUTPUT. 


ASAPl 


(2)  Read  in  namelist  containing  variables  describing  date  and  time  and 
whether  this  is  the  first  analysis  in  a  sequence. 

(3)  Open  unit  2,  the  unpacked  CUE  Level  II  data  set,  described  by  Norquisc 
(1984,  p.  17-33)  and  check  header  record  date/time  against  that  read  via 
namelist  (2) . 

(4)  Set  up  the  vertical  structure  in  /SICK/.  (See  Section  5.4.) 

(5)  Read  In  tvne  1  data  (RAOBS)  from  file  2  on  unit  2  (MASTORl).  In  (5)  - 

(8)  the  call  to  SETFG  initializes  the  horizontal  interpolation 
procedure . 

(6)  Read  in  tvne  4  data  (SATEMS)  from  file  5  on  unit  2  (MAST0R4) . 

(7)  Read  in  tvne  2  data  (AIREPS)  from  file  3  on  unit  2  (MAST0R2). 

(8)  Read  in  type  6a  data  (SATWINDS)  from  file  7  on  unit  2  (MAST0R6) . 

(9)  Count  valid  data  read  in  and  display  totals  (CNTOBS). 

(10)  Gross  and  buddy  check  data  (FLAGS),  assign  points  to  each  observation 
(POINTS) .  These  points  will  be  used  in  the  observation  selection 
procedure.  Count  remaining  data  and  display  totals  (CNTOBS). 

(11)  Store  residuals  on  unit  8  and  identification/location  parameters  on 
unit  3. 

(12)  Beginning  of  analysis  phase  ((12)  -  (27)].  Select  tropospheric  value  of 

(2.0E-12  m'^)  given  by  Dey  and  Morone  (1985).  (See  (23).)  Begin 
loop  over  all  latitudes  in  the  analysis  grid. 

(13)  Determine  ^  ,  the  latitude  of  the  current  grid  row.  The  sine  of  ^  is 

O  O 

read  from  unit  7  if  Guassian  latitudes  are  used  (ITLAT  -  1),  otherwise 
regularly  spaced  ^  are  used.  (Currently  ITLAT  is  hard-wired  to  1 . ) 

O 
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(14)  If  this  is  the  first  analysis  in  a  sequence  read  in  the  standard  esti¬ 
mated  prediction  errors  (EPE)  at  the  mandatory  levels  from  unit  10.  In 

this  case  the  EPE,  E  (2.),  E  (0,),  E  do  not  depend  on  longitude. 

pi  pi  pi 

(15)  Begin  loop  over  all  longitudes.  Define  A 

O 

(16)  If  I  >  70,  convert  (A  ^  )  to  (x  y  )  the  coordinates  in  a  polar 

O  O  O  D  O 

stereographic  projection  true  at  the  pole.  (See  Section  5.5.) 

(17)  Define  the  latitude/longitude  observation  selection  box.  The  box  is 
±  15*  in  latitude  and  ±  30*  in  longitude  if  |^  |  <  60*.  (Longitude 

.o  * 

window  would  be  ±  40*  if  <  1.8  E  -  12  m  ,  but  this  will  not  occur 
for  current  parameter  settings.)  For  \4>g\  >  ^0*  longitudes  will  be 

O 

considered.  (See  (20)  below.) 

(18)  Begin  data  selection  phase  ((18)  -  (25)].  (See  Norquist  1986b, 

pp.  9-10.)  Set  to  0.1.  Start  loop  on  data  type.  NTY  -  1,  2,  3,  4, 

which  correspond  to  type  -  1,  4,  2,  6a.  If  (i.e.  NSMDS) 

observations  have  been  selected  exit  loop.  If  there  are  no  observations 
for  this  type  go  to  end  of  loop  on  type  (25).  Begin  loop  on 
observations  within  this  type .  This  algorithm  implies  that  if  type 

1  observations  are  found  with  other  data  types  are  not 

considered,  but  all  the  type  1  observations  are  considered.  Similarly 
for  type  1  and  4  combined,  type  1,  4,  and  2,  combined,  etc. 

(19)  Check  that  is  within  the  selection  box.  If  not  make  use  of  fact  chat 
data  are  sorted  by  latitude.  If  is  south  of  box  go  to  end  of 
observation  loop.  If  is  north  of  box  go  to  end  of  type  loop. 

(20)  If  l^gl  <  60,  check  that  A^^  is  within  the  selection  box.  If  not  go  to 

end  of  observation  loop.  Note  the  check  when  the  box  straddles 

Greenwich,  i.e.  when  A„,„  <  A„.  .  In  this  case  A„„„  <  A_  <  A„.  implies 
’  max  min  max  n  min 

Aj^  is  not  in  the  box. 
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(21)  If  l^gl  S  70  then  transform  observation  location  (ij^.  to  polar 

stereographic  projection  coordinates  and  calculate  d^^^  the 

distance  between  the  grid  point  and  the  observation. 


(22)  Else  if  |^  |  <  70  then  calculate  distance  with  approximate  great  circle 

D 

formula  (Schlatter,  1975). 


(cos  («i)  AA)^  +  (A^)^ 


z  z  2 

(23)  Calculate  horizontal  correlation  u  -  exp  (-k,d  ).  If  the  number  of 

gn  n  gn 

observations  selected,  M  <  skip  observations  with  /i  <  If 

M  i  Mnax'  observations  with  /i*  <  A**iiiin-  Here  *  indicates  that 

has  been  multiplied  by  the  number  of  points  assigned  to  the  observation 
in  (10). 


(2A)  Insert  the  current  observation  in  the  list  of  observations  for  this 
analysis  point  stored  in  descending  p*  sort  order.  If  M  was  already 
Mjiax*  ^rop  the  last  one  In  the  list  and  reset  ^ax  ' 

set  to  the  last  p*  in  the  list.  The  list  is  maintained  in  two 

arrays;  NAKEY(m)  -  n,  which  points  to  the  observation  arrays  and  MOTTJCi' 
(m)  - 

(25)  End  of  data  selection  loops  on  observation  and  type  (18). 

(26)  Apply  statistical  interpolation  procedure.  (ASAP2) . 

(27)  End  of  loops  over  analysis  longitudes  (15)  and  latitudes  (12). 

(28)  End  of  main  module. 
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Name:  ASAP2  (subroutine) 

Purpose:  Performs  update  of  all  variables  (2,0, t/)  at  all  sigma  layers 

for  a  single  (A  ,4>  )  grid  location.  The  procedure  follows 
S  & 

that  of  Dey  and  Morone  (1985)  and  Bergman  (1979)  as  described 
by  Norquist  (1986b,  pp.  10-18). 

In  brief,  ASAP2  obtains  (a)  the  first  guess  values  and 

(b)  the  expected  prediction  errors  at  the  analysis  point. 
Then,  for  each  layer,  including  a  subsurface  layer,  ASAP2 

(c)  selects  observations  highly  correlated  with  the  grid 
point  variables,  (d)  calculates  the  observation  location-grid 
point  prediction  error  correlations  for  the  selected 
observations  (i.e.,  the  r.h.s.  of  the  normal  equation), 

(e)  calculates  the  matrix  of  observation  location-observation 
location  prediction  and  observational  error  correlations  (the 
same  matrix  is  used  for  2,  0  and  ^) ,  (f)  solves  the  normal 
equations  for  the  nondimens ional  analysis  weights  and 
(g)  calculates  the  corrections  and  expected  analysis  errors. 
Finally,  (h)  z  corrections  are  converted  to  T  corrections, 
and  the  first  guess  is  updated  and  stored  along  with  the 
corrections  and  estimated  analysis  errors. 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  ASAPl 

References:  CHLSKY,  ESOBER,  LOWTMP,  RECALC 

Commons  used:  ASIAO,  AS1AS2,  CCONST,  DCONST,  RESID,  SICK 

Arguments:  IGP  index 

JGP  index 

IMON  month  of  the  analysis 
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I/O  units:  U,  6.  10,  11.  OUTPUT 

Description;  Subscript  and  pointer  usage  in  ASAP2  are  more  complicated 

chan  elsewhere: 


rmbolic  name 


IN.  JN 

NVI,  NVJ 


NN 


NI,  NJ 

LI,  U 

ZB(NI,LI),  etc. 


PH(NI,LI),  etc. 


E 


Meaning _ 

i,  j  index  piece  of  data  which  is  used  in  analysis 

D^(i)  data  type  (1,2,3)  -  (z,u,v)  or  (1,2)  -  (z,(u,v)) 

m(i)  intermediate  pointer  to  n  and  i  once  D^(i)  is  known 

n(l)  -  n(m(i) ,Dj(i) )  observation  index  within  local  array 

i(i)  -  i(m(i)  ,D.j.(i))  layer/level  index 

x(i)  -  x(i(i)  ,n(i)  ,Dj(i) )  is  the  value  of  the  i^  datun^^. 

It  is  stored  in  one  of  three  arrays  for  z,  u  and  v 

separately,  depending  on  the  value  of  D.p(i). 

p(i)  -  p(i(i  ) ,n(i) ,Dj(i))  is  the  value  of  pressure  asso¬ 

ciated  with  the  ith  datum.  It  is  stored  in  one  of  two  arrays 
for  z  and  (u,v)  separately,  depending  on  the  value  of  Dj(i). 

expected  standard  deviation,  e.g., 

Ep(Zj^)  EPE  of  the  first  guess 

Ea^^kg^  EAE  (becomes  E^q(2Lj^)) 


EAE  of  previous  forecast 


ASAP2 


(a)  -  (1)  -  (2)  Obtain  the  first  guess  values 

(1)  From  unit  -  4,  read  in  first  guess  values  (Tj^,  0^^,  P^)  for  current 

location.  Compute  P^^  and  Here  and  below  the  subscript  u  denotes  the 

underground  layer. 

(2)  Calculated  G  ,  the  coefficient  of  geostrophy  at  the  grid  point. 

O 

(See  Norquist  1986b,  p.  14  and  Section  5.6.) 

(b)  -  (3)  -  (7)  Determine  the  expected  prediction  errors  (EPE). 

(3)  Begin  determination  of  EPE  (expected  prediction  errors).  If  this  is  the 
first  analysis  in  a  sequci.-e  (IFRST  -  1),  then  use  the  a  priori  NMC 
values  interpolated  linearly  in  fn  p  to  the  sigma  layer  pressures 
Above  SO  mb,  the  50  mb  value  is  used,  below  1000  mb  extrapolation  is 
used.  Since  the  NMC  values  do  not  depend  on  X  they  are  read  outside  the 

O 

longitude  loop  in  ASAP1(14)  and  stored  in  /AS1AS2/.  The  EPE  values 
Ep(2k) ,  Ep(Ok),  Ep(Vk)  are  stored  in  variable  EPE.  Note  that  Ep(Zj^)  is 
the  EPE  of  Z  at  (7)). 

(4)  Otherwise,  if  this  is  not  the  first  analysis  in  a  sequence  (IFR*^!  -  0), 
then: 

From  unit  -  10,  read  expected  analysis  error  of  the  previous 

analysis.  Interpolate  E^jC^k^  growth  during  the  interval  between 

analyses  from  the  mandatory  layer  values  of  Dey  and  Morone  (1985, 

Table  3,  p.  309).  Also  interpolate  E^(2k).  the  climatological  maximum 
analysis  height  errors  from  Table  2  of  Dey  and  Morone,  and  E^(Uk) , 
the  tropical  wind  EPE  from  Table  4  of  Dey  and  Morone.  In  all  three  cases 
the  interpolation  procedure  is  identical  to  that  used  in  (3).  Next, 
force  EgQ  to  be  in  the  interval  (E^,p,E^).  Then  set 

-  E^g(2k)  calculate  Ep(Oj^)  -  using  equations 

(2.12)  and  (B.ll)  of  Dey  and  Morone  (1985)  as  necessary 


ASAP2 


where  kj^  is  set  as  in  FLAGS(7).  The  blended  values  Eg  are  used  only  if 
l^el  <  25.  (Note  that  typographical  errors  in  (B.ll)  have  been  corrected 

O 

here . ) 

(5)  Initialize  variables,  setting  all  corrections  (z,  .  u,  ,  v,  )  to  7'^ro  and 

Kg’  kg 

all  estimated  analysis  errors  (EAE)  to  the  EPE,  e.g.,  “  Ep(Zj^)  . 

(6)  Loop  over  observations  previously  selected  in  ASAPl  (24),  i.e.,  consider 
observations  in  the  list  NAKEY(i),  i-l,NGBOX.  Extract  data  (pj^,  uj^,  vj^, 
Q(u^,  t  P^.  ZA  and  q(Z^))  from  /RESID/  for  each  of  these  observations 
for  each  k  for  which  data  is  present  and  store  in  local  arrays.  Observa¬ 
tions  with  no  data  present  are  skipped.  Since  layers/levels  with  no  data 
are  skipped,  the  pressures  associated  with  the  data  are  stored.  That  is, 
Ujjjj  is  the  u  residual  at  the  fth  layer  which  had  wind  ^  height  data 
present  in  the  nth  observation,  and  the  pressure  of  U|^  is  p^^-  The 
number  of  layer/levels  extracted  is  stored  in  1^.  Local  arrays  of  dsi^. 
^n’  ^n  also  formed.  If  |^  |  >  70,  (u,v)  are  converted  to  polar 
stereographic  corrdinates.  If  N,  the  number  of  observations  extracted, 
is  zero,  go  to  (49). 

(7)  Assign  the  EPE  for  the  residuals,  E  (zA  ),  E  (u,  )  and  E  (v  )  by  finding 

p  fn  p  In  p  in 

the  grid  point  sigma  layer  pressure  closest  to  pA  for  z  or  p^j.^  for 

m  f  n 

(u,v)  and  assigning  the  corresponding  grid  point  EPE  to  the  observation 

EPE.  e.g.,  E  (u^„)  -  E  (Oj^  ). 

*  *111 

(8)  Initialize  IWFLG.  Begin  main  loop  (8*46)  over  layers,  k-l,K+l.  Layer 
K+1  is  the  underground  layer.  Select  k^^  as  in  FLAGS (7 ) .  is  the 
pressure  at  the  grid  point  (i.e.,  P^^,  or  P^  if  k-K+1). 

(c)  -  (9) -(14)  Selects  the  data  to  be  used  at  this  layer. 

(9)  Loop  over  all  observations  n  -  1,N.  Set  L  -  L^. 

(10)  Calculate  the  distance  from  the  observation  to  the  grid  point. 

For  I I  2:  70”,  the  polar  stereographic  projection  is  used.  Calculate 

G  ,  the  coefficient  of  geostrophy  at  the  observation,  and  ,  the  zz 

& 

horizontal  correlation.  (See  Section  5.6  and  Dey  and  Morone  (1985, 

Eq.  2.8).) 


I 
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(11)  Begin  loop  over  all  layer/levels  in  the  observation  i  -  l.L.  If  a  z 

residual  is  present,  then  add  n  and  i  to  list  of  potential  z  observation 

and  levels  (variables  NZT(ni) ,  LZT(iii),  m-l.NTOTZ).  Calculate  the  vertical 

correlation  between  pA  and  P  using  Eq.  (2.9)  of  Dey  and  Morone 

in  g 

(1985)  .  Obtain  the  estimated  observation  error  (EOE)  for  this 

residual,  E  (za  )  (ESOBER) .  Define  «  as  E  (zA  )/E  (zA  ). 

°zz^^  zu  m  o  in  p  in 

Calculate  p'  ,  p'  and  pL^ ,  the  adjusted  total  correlations  between 
'^mg  '^mg  '^mg 

the  height  residual  and  the  grid  point  variables.  (Note  that  m  implies  a 
particular  i  and  n,  while  g  implies  k.)  These  adjusted  correlations  are 
the  correlations  obtained  following  Oey  and  Morone  (1985,  Eqs.  (2.7), 
(A2),  (A3),  (B2),  (B3)  as  appropriate),  divided  by  (1  +  ■  These 

adjusted  correlations  are  the  normalized  weights  one  would  obtain  for 
updates  made  using  each  residual  separately.  The  sum  of  the  absolute 
values  of  these  weights  is 

-  I-.'"!  +  IP'"'!  ♦  IP'^I 

ra  '^mg  '  ■'^mg  '  ''^mg 

For  the  underground  layer  (i.e.  if  k  -  K+1)  the  exact  same  selection 
process  is  used.  This  ensures  that  the  lowest  layer  and  the  underground 
layer  analyses  are  consistent. 

(12)  Repeat  (11)  for  wind  data.  If  wind  data  are  present,  then: 

Add  n  and  1  to  list  of  potential  (u,v)  observations  and  layers  (variables 
NU(m),  LW(m) ,  m-l,NTOTW).  Calculate  vertical  correlation 
between  p^^^  and  Obtain  define  as  (“fn>- 

(Note  that  S<“in>  is  used  for  v  as  well.)  Calculate 

^mg^’  ^mg^'  ^mg'^’  adjusted  correlations,  as  in 

(11),  but  here  using  Eqs  (A4)-(A9)  or  (B^^)-(B9)  of  Dey  and  Morone 
(1985).  Set 

(2)  ,  ,uz,  ,  ,vz,  ,  ,uu,  ,  ,uv,  ,  ,vu,  ,  ,vv, 

m  mg  '  mg  '  mg  mg  '  '  mg  '  mg  ' 

End  of  loop  (9). 


(13)  Determine  the  NLRHS  (currently  10)  largest  )pl  greater  than  0.1  from  the 

combined  set  of  p^^^  and  p^^^  .  At  this  point,  a  (u,v)  wind  is  considered 
ffl  m 

a  single  piece  of  data.  Pointers  to  the  selected  residuals 
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are  saved  in  m(i),  D^(zi),  i-1 , I  (variables  NNRllS(i),  NVRHS(i), 
i— 1,NKT).  The  insert  sort  algorithm  (see  ASAP1(2^))  puts  these  pointers 
in  the  list  in  descending  |/>|  order. 

(lA)  If  I-O  and  k  <  K,  go  to  (A5).  Else  if  I-O  and  k-K+1,  go  to  (A6) . 

Else  I>0;  continue  with; 


(d)  -  (15)  -  (20)  Calculate  the  observation- gridpoint  prediction  error 
correlations  for  the  selected  residuals.  At  this  point,  the  u  and  v  wind 
components  begin  to  be  treated  separately.  Note  that  in  (17)- (19)  the 
correlations  are  not  adjusted  by  dividing  by  (1  +  <„  )  ss  in  (11)  and  (12). 
Note  that  all  correlations  for  the  underground  layer  are  calculated. 


(15)  Set  j“0.  Loop  over  selected  data  i-1,1  Extract  D^(i)  and  n(i)  and  i(i) 
from  m(i) . 

(16)  Repeat  (10),  calculating  d  ,  G  and 

—  ng  n  mg 

(17)  If  D^(i)  -  1,  then  the  current  residual  is  a  height  residual:  Increment 
j.  Add  it  to  the  list  of  pointers,  by  setting  m(j)  -  m(i)  and  Dj(j)  -  1 

22  211  ZV 

(variables  NNRHSV,  NVRHSV) .  Repeat  (11)  but  store  p  ,  p  ,  and  p  ^ 

ng  ng  ng 

as  the  three  r.h.s.  of  the  normal  equations  (variable  (RHSV) . 

(18)  If  Dj(i)-2,  then  the  current  residual  is  a  wind  residual: 

Add  the  u  component  to  the  list  of  pointers  as  in  (17) .  Calculate 

and  as  in  (12).  Store  these  adjusted  correlations  as  the 
mg  mg  mg 

three  r.h.s.  of  the  normal  equations. 

(19)  If  D<p(l)-2,  repeat  (18)  for  the  v  component.  Note:  Dj(j)  -  3.  This  is 

the  point  where  a  wind  residual  is  treated  as  two  pieces  of  data, 
vz  vu  w 

Calculate  p  ,  p  ,  and  p  as  in  (12)  but  store  as  the  three  r.h.s.  of 
mg  mg  mg 

the  normal  equations. 

(20)  End  of  Loop  (15).  Set  I-j .  The  r.h.s.  and  the  pointers  Dj(i),  m(i)  are 
complete . 


(e)  -  (21) -(32)  Calculate  the  matrix  of  observation-observation  correlations. 


(21)  Loop  on  all  pieces  of  data  i-1,1.  Extract  D.j.(i),  n(i),  f(i),  p(i)  and 
Ep(x(i)).  Obtain  Eg(x(i))  (ESOBER)  and  set  -  Eq(x( i ) )/Ep(x( i) ) . 
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(22)  Calculate  variables  needed  for  the  distance  calculation  and  the 
correlation  calculation  which  depend  on  the  i^  piece  of  data.  Calculate 
Gj ,  the  coefficient  of  geostrophy. 

(23)  Loop  on  all  pieces  of  data  j-i.I.  If  j-i,  go  to  (31).  Otherwise,  o'-tain 
D^-Cj),  n(j),  l(j).  p(j).  Ep(x(j)).  E^(x(j))  and  fj  as  in  (21). 

(24)  Calculate  the  distance  d^j ,  the  coefficient  of  geostrophy  Gj ,  and  other 
variables  needed  for  the  correlation  calculation. 

(25)  Calculate  variables  for  the  vertical  correlation  calculation.  Calculate 

zz 

the  horizontal  correlation  u. . . 

U 

(26)  Set  horizontal  observation  error  correlation  (Norquist,  1986b,  pp.  14-16) 
as 


1 


exp(-k^0‘*^^ 


if  sane  variable  and  same  obseirvation 

“  n(i)) 

for  type  1  or  2 
for  satellite  heights 


0 


K  Otherwise 


Note  that  satellite  cloud  drift  winds  errors  are  assumed  to  be 
uncorrelated  horizontally.  In  the  sorted  Level  II  data,  each  aircraft 
report  is  a  separate  "profile"  and  therefore  assumed  to  be  uncorrelated 
here . 


(27)  Calculate  vertical  error  correlations.  Let  x 
Define 


2n^ 


P(i)/P(j) 


Oij 


(1  +  kpQx)-^ 

Interpolate  linearly 
in  X  in  look-up  table 


for  heights,  type  1  or  2 

for  winds,  type  1  or  2 

for  satellite  heights 
(Norquist,  1986b,  Table  3) 


0 


otherwise 
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(28) 


(29) 

(30) 

(31) 

(32) 

(f)  - 

(32a) 


(32b) 


(33) 

(34) 

(35) 


and 

i/.  .  -  (1  +  k  x) 
P 


For  1^  I  >  70,  calculate  p  and  t)  using  formulas  from  appendix  A  of  Dey 
& 

and  Morone  (1985).  Each  of  the  9  cases  zz ,  zu,  zv,  uz ,  uu,  uv,  vz ,  vu, 

w  are  handled  separately  here.  The  statistical  model  for  p  is  also 

used  for  ij,  only  the  constants  and  k  are  replaced 
w 

by  k^Qt  kpQ  and  '^pQ-  Note  chat  the  use  of  the  geostrophic  model  for 
observation  error  intercorrelations  has  no  effect  since  in  every  case 
that  fiQ^Q  is  non-zero,  the  geostrophic  model  reduces  to  a  factor  of  1. 


For  1^  I  <  70,  calculate  p  and  rj  using  formulas  from  appendix  B  of  Dey 
S 

and  Morone  (1985). 


Set  matrix  entries  R. .  -  R. .  -  p..  +  t.t.n... 

2  ^  ^ 

Set  matrix  entry  R,,  -  1  +  <,  . 

11  1 

End  of  loop  (23)  on  j.  End  of  loop  (21)  on  i. 


Go  to  (32)  . 


(33)  -  (39)  Solve  the  normal  equations 


Geostrophic  decoupling  is  effected  for  grid  points  equatorvard  of  11 
degrees.  This  is  accomplished  by  setting  all  height  wind  correlations 
to  zero.  This  makes  the  analysis  of  height  univariate  and  the  analysis 
of  winds  bivariate.  The  correlation  matrix  becomes  block  diagonal, 
corresponding  to  two  separate  coexisting  linear  systems.  This  approach 
could  be  extended  to  decouple  u  and  v  if  desired.  Note  that  D.p(i) 
equals  1,  2,  or  3  if  the  Ith  datum  is  a  z,  u,  or  v  value  respectively. 

Obtain  least  square  solutions.  First  transpose  matrix  of  right  hand 
sides.  Then  call  CHLSKY  to  obtains  solutions.  CHLSKY  iteratively 
checks  the  solution,  rejecting  highly  correlated  data  as  necessary. 

Begin  loop  on  data  type  Dj(g)  -  1,3. 

Unpack  Ag  and  a'  values  for  this  variable. 

(null) 
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(36)  (null) 

(37)  (null) 

(38)  (null) 

(39)  (null) 

(40)  Convert  a'j^  to  the  dlmensionallzed  weights  using  (Norquist,  1986b, 

Eq.  11) 

-  a'j  Eg/Ep(x(i)) 

where  Eg  -  EPE(k ,  D,j.(g) )  for  k  <  K  and  Eg  -  Ep(Z^)  for  the  underground 
layer . 

(null) 

Calculate  C,  the  correction.  C  -  Ea«x(i). 

i  ^ 

Perform  gross  error  check  of  C  against  ZER  (250  m)  or  UER  (25  m/s) 
(Norquist,  1986b,  p.  17).  If  C  passes  the  check,  go  to  (44). 

Set  C-0,  Ag-l .  If  C  is  a  wind,  set  the  flag  (IWFLG(K))  so  that  both  u 
and  V  corrections  are  set  to  zero  in  (48) . 

If  k-K+1,  set  Z^^g-C;  otherwise,  set  E^(Zkg^'  ^a^“kg^  ^a^'^kg^  equal  to 

times  Ep(2i,),  Ep(aj^)  or  Ep(\)  and  u^^  equal  to  C  if 

D-j(g)  -  1,  2  or  3.  (See  Eq.  2.11  of  Bergman,  1979.)  End  of  loop  (33) 
on  Dj(g) .  Go  to  (46) . 

(45)  Set  Zj^g  -  vij^g  -  Vj^g  -  0  and  associated  EAE  to  the  corresponding  EPE. 
(Like (5)  but  for  single  k  only.) 

End  of  main  loop  (8)  on  layers. 

Convert  zj^g  to  tj^g  (LOWTMP)  and  update  first  guess  ”  "Tyt  +  tj^g- 

Update  (0j^t^\^)  by  adding  (*^RgtVjjg);  first  converting  analysis  winds  back 

to  (A,^)  coordinates  if  ^  |  >  70. 

& 


(^1) 

(42) 

(43) 

(43A) 

(44) 


(46) 

(47) 

(48) 

(49) 


Output  Tj^,  Oj^,  to  unit  6  and  E^(Zj^^),  along  with 

and  V, 


kg  '  “a'"kg'’  “a' 'kg' 

^kg‘  '^kg  Note  that  only  E^(zj^g)  will  be  used  for  the 

next  analysis  (see  (4)).  Return. 
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Name : 

Purpose: 

Author; 

Documentation 

Referenced  by 

References ; 

Commons  used: 

Arguments ; 


I/O  units; 

Description: 
Determine  the 
corner  of  the 

y'  -  y-j .  the 

f(x,y) 
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BILINR  (subroutine) 

Interpolates  a  single  field  bilinearly. 
D.  Norquist,  SASC,  1980  -  1986 
R.  Hoffman.  AER.  1986 

FG 

none 

none 


A 

fj^j,  field  array 

(X,Y) 

(x.y)  location  in  grid  units 

NX 

first  A  array  dimension 

NY 

(not  used) 

ANS 

f(x,y).  Interpolated  value 

none 

coordinates  (i.j)  of  the  grid  point  in  the  lower  left  hand 
grid  cell  enclosing  the  point  (x.y).  Determine  x'  -  x-l, 
(x.y)  location  with  respect  to  (l.j).  Then  set: 

‘  *'y'  ^i+i,j+i  +  fi+i,j  +  y'^i.j+i 


+  d-x')  d-y')  f^j 


A-20 


CALCRES 


Name:  CALCRES  (subroutine) 

Purpose;  Calculate  residuals  for  single  level  data. 

Author;  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER.  1986 

Referenced  by:  MAST0R2,  MAST0R6 

References:  None 

Commons  used;  DCONST.  DEL,  FGDATA,  UASIGMA 

Arg\ifflents :  PO 

POH 

I/O  units:  none 

Description: 

(0)  Residuals  are  returned  In  /DEL/, 

(1)  Set  t,  z  residuals  to  missing. 

(2)  Calculate  (u,v)  residuals.  If  (U^  •  )  is  missing  set  (u,v)  residuals 

to  missing.  ^ 
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NAME: 

Purpose : 


Author : 

Documentation : 

Referenced  by: 

References : 

Commons  used: 

Arguments : 

I/O  Units; 


CHLSKY  (subroutine) 

Solves  the  normal  (i.e.  least  squares) 

equations  Ax  -  b  where  A  is  NxN  and  b  is  NxM,  using  the 
Cholesky  R  R  decomposition.  For  the  wind  and  height  analysis 
M  -  3.  This  version  does  no  pivoting  or  shifting  of  data 
arrays.  For  ill-conditioned  cases  highly  correlated 
variables  are  effectively  removed  but  remain  as 
placeholders.  For  example  if  the  ich  variable  is  to  be 
removed  then  the  ith  equation  becomes  x(i)  -  0.  In  addition, 
if  the  maximum  absolute  value  of  x  is  too  large  or  if  the 
computed  value  of  A£  -  1  -  (x,b)  is  not  in  [0,  1]  the  most 
highly  correlated  variable  is  removed. 

Variables  x  and  b  may  share  storage.  Only  the  upper  half  of 
A  is  referenced.  H  must  be  smaller  than  MAXN  and  M  must  be 
smaller  than  MAXM,  so  that  work  arrays  are  large  enough. 

R.  Hoffman.  AER,  1988 

R.  Hoffman,  AER,  1988 

ASAP2 

SPOCO,  SPOSL  (Llnpak  routines) 

None 

A,  b,  N,  M,  NDIH,  x,  Ag  as  described  above.  The  leading 
dimension  of  A,  b,  and  x  are  all  NDIM  in  the  calling  program. 
Note  that  because  of  the  use  of  work  arrays  A  and  b  will  not 
be  changed. 

None 
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Description: 

(0)  Declare  failure  criteria  parameters  and  subroutine  arguments.  Declare 
work  arrays.  Currently  the  maximun  system  size  is  20  and  the  maximum 
number  of  right  hand  sides  is  3.  If  these  sizes  are  exceeded  the  program 
aborts . 

(1)  Initialize  work  arrays  by  copying  A  into  WA  and  b  into  WB. 

(2)  Decompose  A  -  R^R.  Preset  RCOND  to  0.  If  A  is  singular  RCOND  is 
unchanged  by  SPOCO. 

(3)  If  algorithm  failed  because  of  a  small  or  zero  RCOND  go  to  (8).  (Since 
pivots  are  returned  along  the  main  diagonal  of  WA  a  possibly  cheap 
estimate  of  whether  A  is  poorly  conditioned  might  be  obtained  by 
examining  the  pivots.) 

(A)  Solve  for  the  x's  by  backsubstitution. 

(5)  Calculate  the  maximum  value  of  x  and  the  values  of  Ag. 

(6)  If  algorithm  failed  because  of  improbable  value  of  XMAX  or  of  Ag  go  to 

(8). 

(7)  Copy  solutions  in  WB  to  x  and  return. 

(8)  -(10)  Error  handling  section:  In  case  of  an  error  remove  one  of  the 

observations  which  is  most  highly  correlated  with  another.  Then  go  back 
to  start  of  routine.  This  assures  us  that  the  each  solution  makes  use  of 
exactly  the  same  data. 

(8)  Find  IMAX,  the  index  of  the  largest  offdiagonal  element  of  A.  Note  that 
large  negative  correlations  are  unlikely  and  are  not  considered  here. 
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(8.1)  Write  values  to  file  for  later  processing 

(9)  Set  ith  equation  to  Xj  -  0  for  i  -  IMAX. 

(10)  Try  again.  Go  to  (1). 
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Name:  CNTOBS  (subroutine) 

Purpose  Counts  and  prints  out  the  ntunber  of  valid  residuals 

classified  by  level  and  by  type . 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman.  AER,  1986 

Referenced  by:  ASAPl 

References:  none 

Commons  used:  DCONST,  RESID 

Arguments :  NUM  N 

INUM 

ICHK  -  0 

-  1 

I/O  units:  OUTPUT 

Description: 

(1)  Initialize:  write  heading,  zero  counters. 

(2)  Loop  over  all  observations,  switching  type  index  (IT)  as  necessary. 

(3)  At  each  level,  increment  height  counter  if  is  present. 

(4)  At  each  level,  increment  wind  counter  if  (uj^.Vj^)  is  present.  End  of  loop 

(2). 

(5)  Print  results.  Return. 


total  number  of  observations 

index  of  last  observation  of  types  1,  4,  2,  6 

before  consistency  checks 

after  consistency  checks 
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Name;  CQCV  (subroutine) 

Purpose;  Translates  quality  indicators  from  CUE  level  II  data  sets  to 

"generic  quality  marks’  of  Hoffman  (1987)  and  section  2.4  of 
this  report.  (See  Section  5.2.  for  type  1  and  type  2  data.) 

Author:  Ross  Hoffman,  AER,  1988 

Docvunentation:  R.  Hoffman,  AER,  1988 

Referenced  by;  MASTORl,  MAST0R2 

References;  none 

Commons  used ;  none 

Arguments;  IQZ  z,  quality  mark/indlcator 

IQT  T,  quality  mark/indicator 

IQQ  <!•  quality  mark/indicator 

IQW  (u.v),  quality  mark/indicator 

IDSI  -  (not  used) 

ALAT  -  (not  used) 

ALON  -  (not  used) 

P  -  (not  used) 

I/O  units:  none 

Description: 

(0)  On  input  IQZ  contains  the  (3ME  quality  indicator,  on  output  it  contains 
the  ASAP  quality  mark.  The  same  holds  for  IQT,  I()Q  and  IQW.  Establish 
translation  tables  IQH  and  IQV  and  lookup  functions  IC()CV  and  J(X)CV.  The 
maximum  value  of  the  translated  quality  indicators  will  be  used  as  the 
quality  indicator.  Note  IQH(1,  1)  -  0  instead  of  1  to  force  translation 
based  on  the  value  of  IQV  in  this  case.  Similarly  IQV(1,  2)  -  IQV(1, 

3)  “  0.  ICQCV  translates  IH  and  IV,  the  horizontal  and  vertical  quality 
Indicator  stored  in  digits  1  and  2  of  the  input  quality  indicator  and 
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returns  their  maxlmun.  JCQCV  retricts  this  maxiiauin  value  to  the  interval 
[0.991. 

(1)  Determine  type:  standard,  TWOS,  COLEBA  or  UNKNOWN. 

(2)  Translate  each  input  quality  indicator  to  a  quality  mark  in  place  and 
return. 
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Name;  ESOBER  (subroutine) 

Purpose:  Returns  estimated  observation  error  for  a  given  variable  for 

a  given  dsi  at  a  given  pressure  in  a  given  month.  Val'ji''-s  for 
satellite  heights  are  from  Norquist  (1986b,  pp.  36-37);  all 
other  values  are  from  Dey  and  Morone  (1985,  Table  5). 


Author : 

D.  Norquist 

,  SASC,  1980  -  1986 

Documentation: 

R.  Hoffman, 

AER,  1986 

Referenced  by: 

ASAP2 

References : 

none 

Commons  used; 

none 

Arguments : 

IDSI 

dsi 

P 

P,  pressure  of  the  observation 

NV 

variable  index;  1-Z,  2-u,  3-v 

PM 

Pjj,  pressure  at  the  mandatory  levels 

MANDLVL 

L,  the  number  of  mandatory  levels 

EOV 

Eg,  estimated  observation  error 

IMON 

month  of  the  analysis  (not  used) 

Description; 

(1)  For  type  2  data  (AIREPS) ;  for  dsi  -  21,22,  set  E^^  -  U.9  m/s;  for  dsi  - 
23,24,  set  E^j  -  6.0  m/s. 

(2)  For  type  6a  data,  dsi  -  61  (SATWINDS) :  if  p  >  400  mb,  set  E^  -  4.2  m/s; 
if  p  <  400  mb,  set  E^  -  7.5  m/s. 

(3)  For  type  1  data,  dsi  -  13,14,15  (TWOS,  dropwlndsondes) ,  set  E^  -  2.0  m/s. 
For  dsi  -  16  (constant  level  balloon),  set  E^^  -  5.0  m/s. 
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(4)  For  all  others,  interpolate  linearly  in  in  p  in  EOE  table  which  has  four 
categories  (columns)  for  each  mandatory  level  (rows).  If  p  <  50  mb,  the 
30  mb  value  is  used.  If  p  >  1000  mb,  the  1000  mb  value  is  used.  The 
four  categories  are 


Column 

Usape 

dsi 

NV 

1 

RAOB  heights 

11 

1 

2 

RAOB  winds 

11 

2.3 

3 

PIBAL  winds 

12 

2.3 

A 

SAT  heights 

A1 

1 
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Name;  FG  (subroutine) 

Purpose:  Bilinearly  interpolates  first  quess  fields  to  the  observation 

location.  Returns  results  in  /FGDATA/and  in  arguments  2STAR, 
PO,  POH. 

Author;  D.  Norquist,  SASC,  1980  -  1986 

Documentation;  R.  Hoffman,  AER,  1986 
Referenced  by;  MASTORl,  MAST0R2,  MASTR04.  MASTOR6 

References;  BILINR 

Commons  used:  CCONST,  DCONST,  FGDATA,  FGFLDS,  SIGK 

Arguments :  RLAT  ^ 

RLON  A 

ZSTAR  Z* 

PO 

POH  Pa 

k 

IDSI  ^ 

IWAT  flag  for  ft  computation  (see  (3)). 

NRDLAT  number  of  latitudes  which  have  been  read.  (Index 

of  latitude  currently  stored  as  second  latitude  in  /FGFLDS/. 
That  is,  ^2  ■  +  (NRDLAT  -  1)  t^.) 

I/O  units:  1,  5 

Description: 

(0)  Note  that 

(1)  Determine  fine  grid  coordinates  (x,y). 
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(2)  If  y  >  NRDLAT,  values  for  a  new  latitude  are  needed.  Shift  northern 
latitude  values  to  southern  latitude  values  in  /FGFLDS/.  Read  in  a 
new  northern  latitude,  incrementing  NRDLAT.  Go  to  (2). 

(3)  Interpolate  T,  0,  V,  R  from  /FGFLDS/  to  /FGDATA/.  Relative  humidity 
is  interpolated  only  if  data  is  type  1  (RAOBS)  or  if  data  is  type  U 
(SATEMS)  and  I WAT  i  0. 

(4)  Interpolate  P*  from  /FGFLDS/.  Set  and  from  and  /SICK/. 

(5)  If  not  type  6  (i.e.  for  RAOBS,  AIREPS  or  SATEMS)  interpolate  Z*  from 
/FGFLDS/  and  obtain  heights  by  integrating  Tj^  hydrostatically. 
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Name ; 


FLAGS  (subroutine) 


Purpose:  Perform  gross  and  buddy  checks  on  the  data  stored  in  /RES ID/ 

(see  Norquist,  1986b,  pp .  5-8).  All  comparisons  are  made 
within  boxes  approximately  10'  square  (see  Section  5.3). 
Cross  error  check  limit  is  ±3  forecast  error  standard 
deviations.  Buddy  check  follows  Bergman's  (1978,  1979) 
procedure  in  spirit,  with  the  addition  of  keep  flags.  Data 
which  fail  the  checks  are  set  to  missing. 

Author:  D.  Norquist,  SASC,  1980  -  1986 


Documentation;  R.  Hoffman,  AER,  1986 


Referenced  by:  ASAP! 

References:  REJECT 


Commons  used:  AS1AS2,  CCONST,  DCONST,  IFLG,  RESID,  SIGK 


Arguments : 


I  BOX 

list  of  Bj.j 

JDSI 

list  of  dsi 

QLAT 

list  of 

QLON 

list  of 

PTS 

-  not  used 

NWC 

-  not  used 

NUM 

N 

I/O  units:  output 


Description; 

(0)  The  arrays  used  to  store  the  data  flags  and  indices  are  confusing.  In 

the  buddy  check,  within  a  box  at  a  single  layer  or  level,  observations  of 
the  same  type  are  compared  by  pairs  i,  j,  with  i  <  j.  I  will  use  the 
notation  in  this  documentation  that  f!'  means  "i  puts  a  flag  on  j "  ,  since 
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i  is  on  top  of  j  in  the  symbol.  Also  I  use  '^F  for  wind  flags  and  ^F  for 

height  flags.  If  the  data  i  and  j  are  consistent,  then  -  f|  -  0 .  If 

they  are  inconsistent  and  i  is  of  better  quality,  then  -  1  and  F'^  -  0 ; 

if  they  are  of  equal  quality  then  Fj  -  f|  -  1  and  if  j  is  of  better 

quality  then  F-|  -  1  and  f!^  -  0. 

1  J 

In  the  code  the  following  symbols  are  used 


N  -  i 
NO  -  j 


IN  -  i 
INO  -  j-i 


IFWI(IN,INO)  -  ”f| 
IFWJ(1N,IN0)  -  "Fj 


IFZTI  (IN. INO)  -  ^F^ 


IFZTJ  (IN, INO)  -  ^Fj 


(1)  Prediction  error  standard  deviations  ,  Ep('t^) ,  Ep(O^)  and  Ep(V^) 

for  the  mandatory  levels  in  5  latitude  belts  are  defined  from  Table  1  of 
Norquist  (1986b).  The  constants  a  and  b  for  the  buddy  check  (cf.  Eq .  7.1 
of  Bergman  (1979))  are  set  to  3.0  and  1.5,  one  half  of  Bergman's  values. 


(2)  Begin  loop  over  the  18  10”  latitude  bands. 


(3)  Begin  loop  over  the  boxes  in  this  band. 

(4)  Calculate  the  box  index  for  the  current  box  and  begin  loop  over  all 
observations  in  the  box. 


(5)  If  array  overflow  occurs  (current  dimensional  constraint  is  100  observa¬ 
tions)  output  message  and  go  to  (6).  Otherwise,  add  current  observation 
to  the  list  of  observations  in  this  box.  NKB(i)  -  n(i),  i-1,  Ng  is  the 
pointer  to  the  location  in  /RESID/  of  the  i£tl  observation  in  the  current 
box.  End  of  loop  (4) . 

(6)  Loop  over  all  observations  in  the  box  i-1,  Ng,  extracting  residuals  and 
quality  marks  from  /RESID/  and  storing  them  in  local  arrays.  If  |^g|  > 
70”,  winds  are  converted  to  polar  stereographic  winds. 
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(7)  Loop  over  layers  selecting  the  appropriate  kj^  for  this  layer,  (See 

Dey  and  Morone,  1985,  p.  308).  If  S  0.125  kj^  -  2E-12m  if 

0.125  >  Oj^  >  0.070,  then  kj^  -  1.5E-12in'^,  and  if  0.070  >  then 

k^  -  1.0  E-12ni'^. 

(8)  Initialize  all  toss  flags  and  keep  counters  to  zero. 

(9)  Begin  loop  over  all  observations  in  the  current  box.  Determine  indices 

of  the  mandatory  levels  above  and  below  the  layer  and  level  pressures  of 
the  current  data.  That  is  find  ig,  iA  such  that 

Pa  >  ?  a  >  Pa 
I  *  k  i 
b  a 

Pa  &  P  oA  >  Pa 

*  k+1  fA 

Also  determine  the  latitude  belt  for  looking  up  the  forecast  error 
standard  deviations. 

(10)  Interpolate  E  (2a)  to  oa  and  E  (0  ) ,  E  (^  )  to  oi,  linear  in  in  p. 

p  f  k+1  P  P 

Below  1000  mb  extrapolation  is  used,  above  50  mb  the  value  at  50  mb  is 
used. 

(11)  Perform  gross  error  check  on  current  observation  at  current  layer/level. 
Residuals  greater  than  3  forecast  standard  deviations  are  rejected,  by 
changing  their  values  and  associated  quality  marks  to  missing.  REJECT  is 
called  to  report  the  rejected  data  and  the  reason  for  rejection.  A  list 
of  observations  in  the  current  box  which  have  failed  one  or  more  checks 
(gross  or  buddy)  is  maintained;  NCH(m)  for  m-1,  ICH  is  the  mth 
observation  which  has  failed  a  check.  (See  (33)  and  (38).)  Note  that  at 
statement  label  72  NCHT  is  1  if  a  height  error  has  been  detected  and  2  if 
a  wind  error  has  been  detected.  End  of  loop  (9)  on  observations. 

(12)  If  there  is  more  than  one  observation  start  a  loop  over  all  observations 
i-1 ,  Ng.  This  loop  sets  the  initial  flags.  Set  n(i)  -  NKB(i);  this  is 
the  location  of  the  ith  observation  in  /RESID/. 
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(13)  Determine  the  coordinates  of  the  ith  observation  for  distance 
calculation.  A  polar  stereographic  projection  is  used  if  |^g|  >  70. 

(14)  Begin  loop  over  j-  i+1,  Ng.  Set  n( j )  -  NKB(j);  this  is  the  location  of 
jth  observation  in  /RESID/. 

(15)  Calculate  coordinates  of  jth  observation  as  in  (13)  and  ,  the  distance 
between  the  two  observations  as  in  ASAPl  (21-22). 


(16)  Calculate  the  horizontal  zz  correlation  (Dey  and  Morone ,  1985,  Eq .  2.8) 
and  the  vertical  correlation  (ibid.,  Eq.  2.9)  using  layer  pressures. 

If  or  Zj  is  missing  go  to  (19)  .  Otherwise  calculate  the  vertical 
correlation  (Eq.  2.9)  using  level  pressures,  and  the  value  of  of  the 
ZZ  total  correlation  (ibid.,  Eq.  27).  Calculate  the  r.h.s.  and  i.h.s.  of 
the  buddy  check  criterion  for  z  (Bergman,  1979,  Eq.  7.1);  these  are 
(a-bp^j)Ep(x)  and  |x^-Xj|. 

(17)  Determine  Qj^(zj^)  and  Qj^(Zj),  the  z  quality  levels.  (See  Section  5.2.) 


(18)  If  the  buddy  check  criterion  is  violated  for  z,  then  set  toss  flags: 


(a)  If  then  -  1  and 

(b)  If  QL(Zj)  s  QL(Zi)  then  ^fJ  -  1. 


If  the  criterion  is  not  violated  and  if  p^^  i  0.75, 
keep  counters  ^Kj . 


then  increment  the 


(19)  If  either  wind  is  missing  go  to  (22).  Calculate  p^^  and  using 

equations  from  Dey  and  Morone  (1985;  either  (A6)  and  (A9)  or  (B6)  and 
(B9)  as  appropriate) . 


(20)  Determine  wind  quality  levels  ■"''j )  • 

QL(Uj.Vj). 
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(21)  Set  flags  for  winds  as  in  (18).  The  buddy  check  criterion  is  violated  if 

either  u  or  v  components  violate  their  separate  tests.  The  keep  counters 

-  mi  vv 

are  incremented  only  if  p..  and  p. .  are  both  >  0.75. 

ij  iJ 

(22)  End  loop  (14)  on  j.  End  loop  (12)  on  i. 

(23)  Loop  over  all  observations  in  the  box.  Remove  all  toss  flags  on 

observations  which  have  2  or  more  keep  flags,  i.e.  If  (^Kj^  >1)  then 
z  i 

-  0  for  j  -1,  Ng  and  similarly  for  the  winds. 

(24)  Initialize  interative  rejection  cycle.  The  minimum  number  of  toss  flags 
which  require  an  observation  to  be  tossed  (N.^)  is  set  to  2,  unless  there 
are  only  2  observations  in  a  box  in  which  case  it  is  set  to  1. 

(25)  Begin  rejection  cycle.  Initialize  sums  (^S^,  N|() .  Begin  loop  over 

all  observations  in  box,  i-l,Ng.  (Arrays  IQZ  and  IQU  are  used  here  for 
''Sj^  and  ^S^.) 

(26)  Accumulate  the  number  of  toss  flags  on  i  for  j  >  i.  (Process  f|  for 

j  >  i). 

(27)  Accumulate  the  number  of  toss  flags  on  j  by  i  for  j  >  i.  (Process  F^  for 
j  >  i).  Note  (25-27)  are  equivalent  to  creating  the  sums 


(28)  If  <  Nj  and  then  increment  Nj^  the  number  of  observations 

which  are  assured  of  being  kept.  End  loop  (25)  on  i. 


(29)  If  N|^  -  Ng  then  all  observations  remaining  are  acceptable;  go  to  (35). 
Otherwise  there  are  more  to  be  tossed. 


(30)  Determine  minimum  and  maximum  number  of  toss  flags  (min  S^,  max  S^)  and 

the  Index  i^^^  which  corresponds  to  the  observation  with  the  maximum,  for 
both  z  and  (u,v),  separately. 
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(31)  If  max  $£  is  >  A  truncate  to  a  multiple  of  4  (e.g.  7  becomes  A)  for  both 
z  and  <u,v),  separately. 

(32)  Toss  data  (by  setting  the  local  residuals  to  missing)  for  those  i  which 
have  $£  >  the  (truncated)  max  S^.  REJECT  is  called  to  report  the 
rejected  data  and  the  reason  for  rejection.  Maintain  a  list  of  pointers 
to  the  observations  in  the  local  arrays  which  have  been  tossed  during 
this  rejection  cycle.,  e.g.  NZ(m) ,  m-1,  INTZ  contains  the  pointers  to  the 
tossed  z  observations. 

(33)  Update  the  list  of  pointers  to  observations  which  have  failed  one  or  more 
tests  (NCH(m) ,  m-1,  ICH)  to  reflect  changes  made  in  (32)  using  the  NZ  and 
NV  lists. 

(34)  Remove  all  toss  flags  imposed  by  observations  which  have  just  been  tossed 
on  other  observations.  That  is  for  i-1,  Ng  if  >  (truncated)  max  (S^) 
then  set  Fj  -  0  for  J-1,  Ng,  for  z  and  (u,v)  separately. 

(35)  Go  to  (25)  to  start  next  iteration  unless  10  iterations  have  already  been 
performed. 

(36)  Loop  on  observations  in  the  box  i-l,Ng.  Restore  IQZ,  IQW  arrays  from 
/RESID/,  These  arrays  (Q(z) ,Q(u,v))  were  used  to  store  and  ”S£. 

Since  some  observations  have  been  tossed,  set  Q(z)  to  missing  if  the  z 
residual  is  missing,  and  set  Q(u,v)  to  missing  if  the  (u,v)  residual  is 
missing.  Check  if  any  observations  still  might  be  tossed.  If  this 
occurs,  print  message. 

(37)  End  of  loop  (7)  on  layers. 

(38)  Move  changes  back  to  OBD  array,  using  list  NCH.  If  |^g|  >  70  convert 
winds  from  polar  stereographic  back  to  (X,  4>)  coordinates. 

(39)  End  of  loops  (2)  and  (3)  on  buddy  check  boxes.  Print  message  and  return. 
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Name:  LOWTMP  (subroutine) 

Purpose:  Converts  analyzed  residuals  Zj^  to  temperature  residuals  tj^  and 

corrects  the  first  guess  value  (See  Norquist,  1986b, 

p.  18.  pp.  55-58.) 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 
Referenced  by:  ASAP2 

References:  none 

Commons  used:  DCONST,  LOWT 

Arguments:  CORE  Zj^.  layer  analysis 

CORZU  Zy,  underground  layer  analysis 

TG  first  guess  on  input,  updated  analysis  on 

output 

SL  scaled  pressures  pj^ 

SLU  scaled  pressure  p^  in  the  layer  underground 

Note:  all  arguments  are  variables  for  the  grid  point  being 

analyzed. 

I/O  units:  none 

Description: 

(0)  Vertical  indexing  used  in  this  subroutine  is  illustrated  in  the  following 
figure.  The  first  value  is  used  in  (1)  and  (A).  The  second  value  is 
used  otherwise. 
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k  -  K+l.K 


k  -  K,K-1 


cop  of  the  atmosphere 


J 


k  -  3.2 


k  -  2,1 


k  -  1 


surface 


(1)  Compute  T  corrections  between  two  layers  (i.e.,  at  "interfaces"), 

'c  ■  Pk  ■ ^ . 

These  temperatures  must  be  multiplied  by  -g/R  to  become  dimensionalized 
(see  (3)). 

(2)  Reset  index  so  that  tA  -  tA  k  -  1,...,N 

k  k+1  k 

(3)  Apply  Flattery  algorithm,  treating  ca  as  layer  temperatures.  Since  the 

sigma  structure  /SICK/  doesn't  change,  the  required  constants  are  stored 

T  T  -1  T 

in  /LOWT/.  First  form  BTa  ,  then  A  BTa  and  finally  (A  A)  A  BTa. 

Here  Ta  is  the  vector  of  tj\.  Dimens ionalize  the  result  by  -g/R  to  get  T, 
the  layer  temperatures . 

(4)  Set  uppermost  layer  residual  to  zero:  ”  0-  (Note;  layers  include 

subsurface  layer.)  Update  temperatures  ^k  "  ■*'  *^k+l  ^  ~ 
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Name:  MASTORl  (subroutine) 

Purpose;  Reads  tvoe  1  data  (RAOBS)  from  unpacked  CUE  Level  II  data  set, 

interpolates  to  sigma  and  calculates  residuals,  which  are  then 
stored  in  /RESID/.  (See  Norquist,  1986b,  pp.  31-33.) 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation;  R.  Hoffman,  AER,  1986 

Referenced  by:  ASAPl 

References;  CQCV,  FG,  PTOSIG 

Commons  used;  DCONST,  FGDATA,  RESID,  UADATA,  UASIGMA 

Arguments ;  IDATE  date 

ITIME  time 

NW  -  (not  used) 

NUM  N 

QLAT  list  of 

QLON  list  of 

JDSI  list  of  dsi^ 

IBOX  list  of 

I/O  Units:  2, OUTPUT 

Description: 

(1)  Initialize  time  and  date  window  as  ±  3  hours  about  date/time  which  was 
input.  Set  ICHNG  -  0;  this  is  a  flag  used  to  overwrite  a  previous 
observation  with  a  new  observation  of  higher  quality.  Initialize  N,  L. 

(2)  Read  a  record  from  unit  2,  using  format  described  by  Norquist,  (1984, 
pp.  18-19).  On  EOF  go  to  (19). 
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(3)  If  not  acceptable  go  to  (2).  Specifically  dsi  -  75,  cloud  data,  are 
skipped.  The  presence  of  other  non- type  1  data  results  in  an  ABEND. 
Records  with  missing  p,  negative  z  or  date/time  not  in  the  window  are 
skipped . 

(4)  If  X,  or  dsi  do  not  match  or  dsi^j .  go  to  (12).  That  is  if 

current  record  is  not  part  of  current  observation,  it  is  time  to 
complete  processing  the  current  observation  and  begin  new  one.  Note 
that  the  test  makes  use  of  fact  that  A,  ^  are  coded  only  to  hundreths  of 
degree. 

(5)  If  level  dimension  (currently  90)  is  exceeded  in  /UADATA/,  skip  this 
level,  go  to  (2). 

(6)  If  this  is  not  the  first  level  for  the  current  observation,  i.e.,  if 
i  ^  1,  go  to  (8).  Set  Ajj,  dsi^j  to  A,  dsi .  Calculate  (see 
Section  5.3). 

(7)  If  N  >  1,  then  loop  over  all  previous  observations  searching  for  a 

matching  A^.^,  If  match  Is  found  then:  If  dsi  <  dsi^^.  then  the  new 

observation  is  better  (see  WHO  (1978),  Appendix  A,  Table  1).  Set 
ICHNG  -  1  to  force  overwriting  observation  n.  (To  do  this  temporarily 
set  N  -  n.)  Go  to  (8),  Else  if  dsi  s  dsi^.^.  then  flush  current  records 
until  new  location  is  found.  Go  to  (3). 

(8)  Translate  T  in  Celsius  to  T  in  Kelvin,  T^  to  R,  wind  speed  and 
direction  to  (U.V).  If  quality  indicator  is  99  for  any  variable,  set 
corresponding  value  to  missing  (DNN  -  -999.9).  If  T  is  missing,  or  if 
dew  point  or  dew  point  depression  is  unreasonable  set  R  to  missing. 

Uinds  at  the  pole  are  set  to  missing. 

(9)  Combine  quality  indicators  and  translate  into  qualtiy  marks  (CQCV) . 
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(10)  Set  all  quality  marks  >  1  to  99,  set  corresponding  values  to  missing. 
This  keeps  only  data  which  have  been  checked  satisfactorily  and  data 
which  have  not  been  checked. 

(11)  Increment  L.  Go  to  (2). 

(12)  Save  X,  ^  for  the  new  observation  which  has  been  detected.  If  L  -  1 

•  (which  can  happen  only  if  N  -  1) ,  go  to  (6). 

(12a)  Eliminate  duplicate  and/or  out  of  order  levels  from  current  observation. 

(13)  Calculate  first  guess  values  at  for  Z  ,  P  ,  t  ,  0  ,  V 

^  Ic  k  Ic  Ic 

and  2a  (FG) .  Although  IWAT  -  0,  is  calculated,  since  this  is  type  1 
data. 

(14)  If  the  observation  has  only  a  single  level  skip  it  by  going  to  (17a) . 
Interpolate  observation  to  sigma  coordinates  (PTOSIG) .  If  all  sigma 
coordinate  data  is  missing  go  to  (17a). 

(15)  Compute  residuals.  Make  sure  that  missing  data  results  in  missing 
residuals . 

(16)  Store  residuals  and  associated  quality  marks  in  /RESID/.  (At  this  point 
quality  marks  are  either  0  for  checked  residuals,  1  for  unchecked 
residuals,  or  99  for  missing  residuals.) 

(17)  Restore  N,  ICHNG  if  overwriting.  Increment  N  for  next  observation. 

(17a)  If  EOF  go  to  (19).  Otherwise  more  contents  of  current  record  to  start 
of  /UADATA/  arrays,  since  this  data  is  start  of  the  new  current 
observation.  This  resets  L  -  1.  Go  to  (6). 

(18)  Position  unit  2  at  start  of  next  file.  Go  to  (12a). 

(19)  Undo  (17);  N  -  N-1.  Return. 
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Name:  MAST0R2  (subroutine) 

Purpose:  Reads  type  2  data  (AIREPS,  etc.)  from  unpacked  GWE  Level  II 

data  set.  From  each  report  (u,v)  residuals  are  calculated  and 
assigned  to  the  analysis  layer  closest  to  the  reported 
height.  (See  Norquist,  1986b,  pp.  33-34.) 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  ASAPl 

References:  CALCRES,  CQCV,  FG 

Commons  used:  DCONST,  DEL,  FGDATA,  RESID,  UADATA,  UASIGMA 

Arguments :  IDATE  date 

ITIME  time 

NW  -  (not  used) 

NUM  N 

QLAT  list  of 

QLON  list  of 

JDSI  list  of  dsij^ 

IBOX  list  of 

I/O  Units:  2 

Description: 

(1)  Initialize  date/time  window  as  ±  3  hours  about  input  date/time. 

Initialize  N  for  next  observation.  Initialize  /UASIGMA/  humidities 
to  missing  (Qj^) •  Initialize  all  of  /UADATA/  to  missing  (P^,  T^,  U^, 

V^,  Q^,  Z^) .  Initialize  all  of  /FGDATA/  to  missing  (Zj^,  Uj^,  Vj^,  Rj^^)  • 
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(2)  Read  a  type  2  header  record.  See  Norquist,  (1984,  pp.  20-21).  At  EOF 
go  to  (23). 

(3)  Check  if  data  is  acceptable.  Non- type  2  data  results  in  an  ABEND.  Data 
outside  date/tine  window  are  skipped  by  going  to  (22). 

(4)  Set  X^,  and  dsi^j  to  values  read.  Calculate  (see  Section  5.3). 

(5)  If  current  observation  is  collocated  with  any  previous  observation, 
including  those  of  other  data  types,  the  current  observation  is  skipped 
by  going  to  (22). 

(6)  Store  Z  and  T  and  associated  quality  marks  from  the  header  record, 
converting  degrees  Celsius  to  degrees  Kelvin. 

(7)  Loop  on  wind  records  associated  with  current  header  record  which  was 
read  in  (2).  (Currently  sorted  data  always  have  a  single  wind  per 
header  record.)  First  wind  report  is  preset  to  missing,  then  data  are 
skipped  if  quality  mark  indicates  wind  information  is  missing  or  if 
location  is  at  the  pole.  Otherwise  wind  speed  and  direction  are 
converted  to  (U,  V)  and  the  quality  mark  is  saved. 

(8)  Skip  the  optional  record  if  present. 

(9)  Translate  the  quality  indicators  to  quality  marks  (CQCV) . 

(10)  Only  quality  marks  of  0  and  1  are  retained.  If  Q(Z)  >  1,  skip  this 
observation  by  going  to  (2).  If  Q(T)  >  1  or  Q(U,V)  >  1,  the  associated 
values  and  quality  marks  are  set  to  missing. 

(11)  Interpolate  first  quess  to  (FG) .  Note  that  R  is  not 

interpolated. 

(12)  If  Z  <  Z  or  Z  >  2a,  then  go  to  (21).  Otherwise  determine  the  index  k^ 

^  K 

of  the  2a  closest  to  Z. 
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(13)  Calculate  the  pressures  associated  with  and  Z  using  formulas  for  a 


Q 

standard  atmosphere  as  described  by  Norquist  (1986b,  Eq .  24  and 

A 

subsequent  discussion) .  Call  these  P(Z)  and  P(Z) . 


(14)  The  actual  pressure  associated  with  the  observation  is  then  taken  to  be 

A 

p  defined  so  that  p  -  P^  P(Z)  -  P(Z) . 

c 

(Note  that  zA  is  at  pressure  pA  .) 

k  ^  ‘^k+l 

(15)  Determine  the  sigma  layer  k^  which  has  pressure  closest  to  p. 


(16)  If  T  is  present,  interpolate  first  guess  linear  in  in  p  to  p  and  then 
add  residual  to  to  obtain  This  effectively  assigns  the 

residual  to  layer  k^.  If  T  is  not  present  set  Tj^  to  missing. 


(17)  If  (U,V)  is  present  interpolate  first  guess  linear  in  in  p  to  p  and  then 
add  the  residual  to  (0.  ,V.  )  to  obtain  (U.  ,V.  ).  This  effectively 

■'v  “Sr  *S^  *S/ 

assigns  the  residual  to  layer  k^.  If  (U,V)  is  not  present, 
set  to  missing. 

(18)  For  all  k  ^  Ky,  set  T^^,  Vj^  to  missing. 

(19)  Compute  residuals  (CALCRES) . 


(20)  Store  results  in  /RESID/.  Each  observation  contains  only  P  ,  Uj^  ,  v 


Increment  N. 


K' 


(21)  If  only  a  single  wind  is  associated  with  the  current  header  go  to  (2). 
This  will  always  be  the  case  with  sorted  data.  Otherwise  process  other 
wind  data  roughly  following  steps  (4),  (11),  (12),  (13),  (14),  (15),  (17), 
(19),  (20).  Then  go  to  (2). 

(22)  Skip  wind  record(s)  and  the  optional  record  if  it  is  present,  which  are 
associated  with  current  header  record. 


(23)  On  EOF  position  unit  2  after  the  EOF  mark.  Undo  last  N  increment  and 
return . 
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Name:  MASTORA  (subroutine) 

Purpose:  Reads  in  SATTEIMS  from  GWE  level  II  unpacked  data  set  which  does 

not  include  SATTEMS  over  land.  Data  is  transformed  to  sip^ma 
coordinate  height  residuals  and  stored  in  /RESID/.  (See  Norquist 
1986b,  pp.  3A-36.)  U.S.A.  Special  Effort  quality  controlled 
and/or  corrected  satellite  temperature  profiles  (dsi  -  41)  are 
preferred  to  raw  reports . 


Author : 


Documentation : 


Referenced  by: 

References : 


Commons  used: 


D.  Norquist,  SASC,  1980  -  1986 
R.  Hoffman,  AER,  1986 

ASAPl 

FG,  SATLTMP 

DCONST,  FGDATA,  RESID,  UADATA,  UASIGMA 


Arguments : 


I  DATE 

date 

ITIME 

time 

NU 

- 

(not  used) 

NUM 

N 

QLAT 

list 

of 

^n 

QLON 

list 

of 

^n 

JDSI 

list 

of 

^n 

I  BOX 

list 

of 

Bn 

I/O  Units:  2 


Description : 


(1)  Initialize  date/time  window  as  ±  3  hours  about  date/time  as  input. 

Increment  N  for  next  observation.  Initialize  observations  interpolated  to 
sigma  as  missing  for  U,  V  and  Q.  i.e.,  DNN .  Initialize 
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all  first  guess  values  at  observation  location  to  missing;  i.e., 

Za  -  0  -  V  -  R  -  T  -  DKN. 

k  k  k  k  k 

(2)  Read  in  type  U  data  record.  This  should  be  the  group  header  record  for  a 
number  of  layer  records  (see  Norquist,  1984,  pp.  25-27).  On  EOF  go  to  (18). 

(3)  If  not  acceptable  go  to  rejection  loop.  Specifically  dsi  must  be  41  and  the 
date/time  must  be  acceptable.  Note:  soundings  over  land  are  not  included 
in  the  input  data  set. 

(4)  Set  Ajj,  dsijj  to  values  read.  Calculate  Bjj. 

(5)  If  collocated  with  a  previous  observation,  possibly  of  another  type,  go  to 

rejection  loop  (17).  However  if  colocated  with  an  observation  of  the  same 

type  and  if  the  data  has  been  corrected  by  the  U.S.A.  Special  Effort  (ISCR 

-  5)  or  found  correct  or  probably  correct  by  the  U.S.A.  Special  Effort 
(IQCI  -  6  or  7)  then  the  new  observation  is  processed.  When  processing  is 
complete  instead  of  incrementing  N  we  shift  the  contents  of  the  Nth  (new) 
observation  to  the  colocated  observation  (see  (16)). 

(6)  Translate  quality  indicator  (IQCI)  to  quality  mark  according  to  WMO  (1986), 
Appendix  A,  Table  38).  If  ISCR  (WMO  (1986),  Appendix  A,  Table  37)  indicates 
observation  was  corrected  set  quality  mark  to  2.  All  observations  with 
quality  marks  greater  than  1  are  skipped  by  going  to  rejection  loop  (17). 

(7)  Begin  reading  layer  data  records.  See  Norquist  (1984,  pp.  26-27)  for 

details  on  the  four  types  of  layer  information.  MAST0R4  assumes  that 

thickness  and  mean  temperatures  are  not  both  present  and  that  standard  and 

nonstandard  layer  precipitable  water  are  not  both  present.  Obtain  first 

guess  values  for  Z  ,  P  ,  Pa,  T  ,  0  ,  ^  ,  2a,  and  if  layer  precipitable  water 
*  k  k  k  k  k  k 

is  present,  R  . 

k 
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(8)  If  no  thickness  data  is  present  go  to  (12).  Read  dat<i  iu^o  arrays  of  P| 
and  Zj  where  Pj  is  the  pressure  at  the  top  of  layer  i,  is  the  height 
at  relative  to  Pj^.  Thus  -  0.  The  layer  temperatures  T^  are  then 
calculated  hydrostatically  (see  Section  5.1.2).  Go  to  (10). 

(9)  For  mean  temperature  data  P^  and  T^  are  read  in.  Data  associated  with 

missing  pressures  or  with  layers  of  zero  or  negative  mass  (4p)  are 
skipped.  The  pressure  at  the  lower  boundary  should  match  the  previous  upper 
boundary  pressure.  If  it  is  too  large  an  error  exists,  if  too  small,  then  a 
layer  with  missing  T^  is  defined  between  the  two  pressures.  Indexing  used 
here  is  normal:  T^  is  the  mean  temperature  between  Pj  and  • 

(10)  Interpolate  T.  to  T,  (at  sigma  layers)  and  Za  (at  sigma  levels)  using 

X  K  iC 

Flattery  algorithm  and  hydrostatic  relationship  along  with  some  first  quess 
information  (SATLTMP) .  Note  dimensional  limit  of  20  layers  in  matrix  ATA. 

(11)  Set  quality  marks  for  Za  and  T  equal  to  the  quality  mark  defined  in  (6)  for 

k  k 

all  levels  which  have  data  and  to  99  for  missing  data.  The  layer  by  layer 
quality  marks  of  the  Level  II  data  are  ignored.  Go  to  (13). 

(12)  Thickness  data  is  not  present;  set  all  quality  marks  and  all  T^^,  and 
Z^  data  values,  to  missing  in  /UASIGMA/. 

(13)  Skip  records  corresponding  to  layer  precipitable  water  data.  Set  all 
associated  quality  marks  and  data  values  to  missing  in  /UASIGMA/.  If  mean 
temperature  data  is  to  be  read  next,  then  go  to  (9). 

(14)  Begin  computation  of  residuals.  Store  P*  in  /RESID/.  Code  is  present  here 
to  convert  Q  to  R  (see  Section  5.1.4).  However,  Q  will  always  be  missing 
(see  (13)). 

(15)  Calculate  residuals;  if  data  is  missing  in  /UASIGMA/,  the  quality  marks  and 
residuals  are  set  to  missing.  Wind  quality  marks  and  residuals  are 
automatically  set  to  missing.  All  residuals  are  stored  in  /RESID/. 
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(16)  If  new  point  colocates  with  a  previous  point  move  data  there.  Otherwise 
increment  N.  Go  to  (2). 

(17)  Rejection  loop:  Read  thru  all  records  of  the  current  sounding. 

Go  to  (2).  Unexpected  EOF  results  in  an  ABEND  here. 

(18)  EOF  handling  (see  (2)). 


MAST0R6 


Name:  MAST0R6  (subroutine) 

Purpose:  Reads  in  type  6a  data  (SATUINDS)  from  unpacked  GWE  Level  II  data 

set.  From  each  report,  (u.v)  residuals  are  calculated 
assigned  to  the  sigma  layer  closest  to  the  reported  pressure. 
(See  Norquist,  1986b,  pp .  41-42.) 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 

Referenced  by:  ASAPl 


References:  CALCRES,  FG 

Commons  used:  DCONST,  DEL,  FGDATA,  RESID,  UADATA,  UASIGMA 


Arguments ; 


IDATE 

date 

ITIME 

time 

NW 

-  (not  1 

used) 

NUM 

N 

QLAT 

list  of 

^^n 

QLON 

list  of 

JDSI 

list  of 

dsi„ 

I  BOX 

list  of 

»n 

I/O  Units:  2,  OUTPUT 


Description: 

(1)  Set  up  quality  mark  translation  tables  ICP  and  ICU  (see  (7)).  Initialize 
date/time  window  as  ±  3  hours  about  input  date/time.  Increment  N  foj  next 
observation. 

(2)  Initialize  /UASIGMA/  humidities  to  missing  (Qj^)  Initialize  all  /UADaTA/  to 

missing  (P^,  Tj,  U^j,  V^,  Q^,  Z^).  Initialize  all  /FGDATA/  to 

missing  (Za,  t  0  ,  v  ,  R  ) 
k  k  k  k  k 
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(3)  Read  a  type  6a  record.  See  Norquisc  (1984,  pp.  29-30). 

(4)  Check  if  data  is  acceptable.  Non- type  6a  data  results  in  an  ABEND.  Data 
outside  date/time  window  are  skipped  by  going  to  (3). 

(5)  Set  Ajj,  dsi^j  to  values  read.  Calculate  Bjj  (see  Section  5.3). 

(6)  If  current  observation  is  collocated  with  any  previous  observation, 
including  an  observation  of  another  data  type,  the  current  observation  is 
skipped  by  going  to  (3).  Note  collocation  criterion  agreement  to  0.1 
degrees  in  latitude  and  longitude  since  for  these  data  location  is  coded  to 
0.1  degrees,  but  other  data  are  coded  to  0.01  degrees.  However  if  colocated 
with  an  observation  of  the  same  type  and  if  the  IQC  -  6  or  7  indicating  the 
U.S.A.  Special  Effort  found  the  new  observation  correct  or  probably  correct, 
the  new  observation  is  processed.  The  new  point  overwrites  the  old  point 
(see  (18)). 

(7)  Set  quality  marks  Q(T)  and  Q(U,V)  from  quality  indicators  IQCP  and  IQC  read 
at  (3)  using  the  translation  tables  ICP  and  ICW.  These  tables  translate  WHO 
(1986)  Appendix  A  Tables  30  and  36  into  the  ASAP  "generic  quality  marks" 

(see  Section  5.2). 

(8)  If  T  is  missing  or  Q(T)  >  1,  set  T  and  Q(T)  to  missing.  Otherwise  convert 
degrees  Celsius  to  degrees  Kelvin. 

(9)  If  wind  speed  or  direction  is  missing  or  Q(u,v)  >  1,  set  (u,v)  and 
Q(u,v)  to  missing.  Otherwise  convert  wind  speed  and  direction  to  (u,v). 

(10)  If  T  and  (u,v)  are  missing  skip  this  observation  by  going  to  (3). 

(11)  Interpolate  first  quess  to  Ajj,  (FG).  Note  that  R  is  not  interpolated. 
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(12)  If  reported  pressure  P  is  not  between  top  and  bottom  sigma  layer 

pressures,  (P  )  then  skip  the  current  observation  by  going  to  (3). 

1  k 

Otherwise  determine  Ky,  the  index  of  the  sigma  layer  pressure  closest  to  the 
observation  pressure. 


(13)  If  T  is  present,  interpolate  first  guess  linear  in  In  p  to  P  and  then  add 
residual  to  Tj^  to  obtain  This  effectively  assigns  the  residual  to 

t-  ir  Jr  j  -  _ Y _ i _ _ nr  _ _ _ x _ X _ 


k^.  If  T  is  missing,  set  Tj^  to  missing. 


(14)  If  (u,v)  is  present  interpolate  first  guess  linear  in  in  p  to  P  and  then  add 

the  residual  to  (Oj^  ,  )  to  obtain  (Uj^  ’  '^k  ^  ‘  effectively  assigns 

the  residual  to  k^.'^  If  Tu,v)  is  missingY  set^(Uj^, Vj^)  to  missing. 

(15)  For  all  k  ^  k^,  set  Tj^,  u^^,  to  missing. 


(16)  Compute  residuals  (CALCRES). 


(17)  Store  results  in  /RESID/.  Each  observation  contains  only  P  ,  u  .  v  . 

*  iC  K 

V  V 

(18)  If  new  point  colocates  with  a  previous  point  move  data  there  (see  (6)). 
Otherwise  increment  N  and  go  to  (3)  unless  N  now  exceeds  6000,  the  current 
array  dimension  in  /RESID/. 


(19)  EOF  or  dimensions  exceeded.  Close  unit  2.  Undo  last  N  increment.  Return. 
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Name ; 

POINTS  (subroutine) 

Purpose : 

Calculates  the  number  of  data  quantity  points  (PTS)  for  each 

observation.  Data  quality  of  SATEMS  is  accounted  for. 

Additional  credit  is  given  if  z  and  (u,v)  are  both  present. 

(See  Norquist,  1986b,  pp.  8-9.) 

Author ; 

D.  Norquist,  SASC,  1980  -  1986 

Documentation : 

R.  Hoffman,  AER,  1986 

Referenced  by; 

ASAPl 

References ; 

none 

Commons  used; 

DCONST,  RESID 

Arguments ; 

JDS  I  dsl^ 

PTS  paints^ 

NUM  N 

IMON  month  index 

I/O  units; 

none 

Description; 

(1) 

Loop  over  all  observations  n-1 .  N.  Initialize  points^  -  0. 

(2) 

Set  factor  to  1.0  unless  dsi^/10-^  (i.e.  SATEMS)  in  which  case 

factor  is  set  to  the  ratio  of  the  observational  standard 

deviations  of  RAOBS  to  SATEMS  (0.62  in  winter  and  0.53  in 

summer) . 

(3) 

Loop  over  all  layers  k  -  1,  K. 

If  ZA  is  present,  add  factor  to  noints... 
kn 

If  is  present  add  factor  to  noints^.^. 

If  zA  and  (u.vlj^j.^  are  oresent  add  factor  to  nointSj^. 

(^) 

End  loop  (3)  on  k.  End  loop  (1)  on  n.  Return. 
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Name ; 


PTOSIG  (subroutine) 


Purpose:  Interpolates  RAOBS  to  sigma.  Assumes  that  temperature  varies 

linearly  in  in  p  and  hence  that  height  varies  quadratic;.! iy  in 
in  p.  Both  and  observations  are  used  to  calculate  ,  the 
sigma  layer  temperatures.  The  T^^  supplemented  by  first  guess 
heights  as  needed  are  used  to  calculate  the  Za.  Winds  and  in 
(specific  humidity)  are  interpolated  linearly  in  inp.  (See 
Norquist  1986b,  pp.  32-33,  61-63.) 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER,  1986 


Referenced  by :  MASTORl 

References:  none 


Commons  used:  DCONST,  FGDATA,  UADATA,  UASIGMA 

Arguments :  NLVL  L 

ZSTAR  Z* 

PO 

POH 

I/O  units:  none 

Description: 

(1)  Begin  loop  over  all  sigma  levels  k-1,  K+1.  This  loop  interpolates  Z^  and  T^ 
to  a  preliminary  estimate  of  Z|\.  Find  ij^,  the  first  i-1,  L  such 

that  ^  Go  to  (3). 

(2)  Fall  through  loop.  If  P^^  is  within  0.1  mb  of  Pa  allow  extrapolation  of  z 
by  setting  -  L.  Otherwise  go  to  (8). 
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(3)  If  -  1  go  to  (8).  (No  dovmward  extrapolation.) 

(4)  If  -  P^)  >  300.  go  to  (8). 

(5)  Find  ij,  the  first  level  i  at  or  below  ij^-1  which  has  both  and  data. 
Find  iy,  the  first  level  f  at  or  above  which  has  both  and  data.  If 
either  can't  be  found,  go  to  (8). 

(6)  If  (P,  -  P,  )  >  300  mb,  go  to  (8). 

(7)  Interpolate  Z^  and  quadratically  in  in  p  to  obtain  the  preliminary 
estimate  of  Za.  Note  that  here  Z^  is  the  height  at  not  at 

This  algorithm  is  comprised  of  Eqs.  C7 ,  C3 ,  C6 ,  C2  and  Cl  of  Norquist 
(1986b,  pp.  61-62),  except  that  the  constants  A  and  B  are  multipled  by  (-1) 
everywhere.  The  quality  mark  Q(Z^)  is  set  equal  to  the  maximum  of  the 
quality  marks  associated  with  the  four  input  data, 

T^  and  Z^  at  and  i^.  Go  to  (9). 

(8)  Do  not  interpolate.  Set  Z|S  to  missing. 

(9)  End  of  loop  (1)  on  levels. 


(10) 


Compute  Tj^,  k-l,K  hydrostatically  from  preliminary  estimates 

of  Z|\  calculated  in  loop  (l)-(9).  Set  Q(Tj^)  equal  to  the  maximum 

of  Q(Za)  and  Q(Za  ).  If  either  of  Za  or  Za  is  missing,  then  T 
k  k^  1  k  k^  X 

missing. 


k 


is 


(11)  Calculate  final  for  k-l,K  hydrostatically  from  Tj^  and  Z^.  Now 

Z|\  is  at  Zjj  is  an  estimate  of  Z  at  the  bottom  of  layer  k.  For 

k-1,  Zjj  -  Z  ;  for  k>l,  Z^-  ^  i  present  or  -  Za  if  Za 

^  k  *1  k”l  k*l 

missing.  If  Tj^  is  missing,  then  Z^  is  missing.  Quality  marks  are  set 

that  Q(Z|\)  -  Q(Tj^)- 


is 

so 


(12)  Begin  loop  on  layers  k-1,  K.  This  loop  interpolates  (u,v)  and  fn  q.  Find 
2^,  the  first  £-1,  L  such  that  P^  < 


Go  to  (15) . 


PTOSIG 


(13)  Fall  through  loop.  If  Pj^  is  within  0.1  mb  of  Pj^ .  use  the  values  and  quality 

marks  of  Pj^,  for  V^,  Otherwise, 

(14)  Do  not  interpolate.  Set  missing.  Go  to  (21). 

(15)  If  f.-l  go  to  (14).  If  (P,  ,  -  P»  )  >  300  mb,  go  to  (14). 

1  Zf  1.  /^i 

(16)  Find  i  ,  the  first  level  at  or  below  -f^-l  which  has  (U,V)  data.  Find  fu  the 
first  level  at  or  above  which  has  (U,V)  data.  If  either  can't  be  found 
go  to  (19). 

(17)  If  (P,  -  P,  )  >  300  mb,  go  to  (19). 

f  u 

(18)  Interpolate  winds  at  and  linearly  in  fn  p  to  to  obtain 

Set  Q^^k’'^k^  maximum  of  ^tid  Go  to  (20). 

(19)  Do  not  interpolate.  Set  to  missing. 

(20)  Repeat  (16)-(19)  for  relative  humidity.  (See  Mitchell,  1985.) 

(21)  End  of  loop  (12)  on  layers.  Return. 


REJECT 


Purpose ; 


REJECT  (subroutine) 


This  subroutine  writes  a  record  to  the  REJECT  file 


(unit  9) . 


Author : 


R.  Hoffman,  AER,  1987 


Documentation:  R.  Hoffman,  AER,  1988 


Referenced  by:  FLAGS 


References:  None 


Commons  used:  AS140,  DCONST,  RESID,  SIGK 


Arguments:  IFLAG  Rejection  indicator. 

K  k,  vertical  level  of  rejected  data. 

N  n,  observation  index  of  rejected  data. 

The  input  parameter  IFLAG  is  a  two  digit  integer.  The  first 
digit  indicates  if  the  rejection  is  caused  by  the  height  (1)  or 
wind  (2)  data.  The  second  digit  is  indicates  if  the  rejection  is 
caused  by  the  buddy  check  (1)  or  the  gross  check  (2).  E-g-  21 
indicates  a  (u,v)  datum  rejected  by  the  buddy  check,  k  and  n 
point  to  the  location  of  the  rejected  data  in  the  OBD  array. 

I/O  Units;  9 


Description: 


(0)  When  first  called  open  REJECT  file. 


(1)  Determine  pressure  associated  with  rejected  data. 


(2)  For  the  nth  observation  write  6,  dsi .  IFLAG.  K,  pA  ,  ZA  , 

-  ‘^k+1  k 

Q(za)  ,  p^,  Uj^,  Vj^,  Q(U|^,  Vj^)  to  the  REJECT  file  on  unit  9.  Return. 
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Name;  SATLTMP  (subroutine) 

Purpose;  Computes  sigma  layer  ten.peratures  and  sigma  level  heights  from 

satellite  observed  layer  temperatures. 


Author : 


Ross  Hoffman,  AER,  1988 


Documentation;  R.  Hoffman,  AER,  1988 


Referenced  by:  MAST0R4 

References;  ANCHOR 


Commons  used;  DCONST,  FGDATA 


Arguments ; 


I/O  Units: 


p 

Pa,  level 

T 

Tj0,  layer 

ATA 

workspace 

NLVL 

L 

ZS 

Za 

k 

TSIG 

Tk 

PSIG 

P 

1- 

POH 

.  K 

ZSTAR 

z* 

OUTPUT 

pressures 

temperatures 


Description: 


(0)  Print  message  Identifying  version  when  first  called. 

(1)  Fill  output  arrays  with  missing  values.  Calculate  the  log  of  the  pressures. 

(2)  Determine  the  first  sec  of  contiguous  non-missing  values  (from  K  to  L-1). 

(3)  Call  ANCHOR  to  calculate  T.  and  Za. 

k 

(4)  Print  diagnostic  output  when  first  called.  Return. 
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Name:  SETFG  (subroutine) 

Purpose:  Initializes  the  bilinear  interpolation  by  reading  in 

southernmost  two  latitudes  of  the  first  quess  grid, 
are  sorted  by  latitude  only  two  latitudes  are  needed 
time . 

Author:  D.  Norquist,  SASC,  1980  -  1986 

Documentation:  R.  Hoffman,  AER ,  1986 

Referenced  by:  A3AP1 
References:  none 

Commons  used:  CCONST,  DCONST,  FGFLDS 

Arguments :  none 

I/O  units:  1,  5 

Description: 

(1)  Rewind  units  1,  5. 

(2)  For  latitudes  1  and  2.  read  T,  U.  V,  ft,  for  all  longitudes  and  layers 
from  unit  1.  Similarly,  read  Z*  for  all  longitudes  from  unit  5. 

Note:  File  structure  is  A,  k,  while  /FGFLDS/  stores  data  by  A,  k. 


t  ho 

Since  data 
at  any  one 
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^ .  Input/Output  Usage . 

Program  ASAPl  uses  units  1-8,  10,  11,  INPUT  and  OUTPUT.  Table  4,1 
describes  these  files  briefly.  Table  4.2  lists  the  I/O  activity  in  ASAPl. 

Table  4 . 1 

Files  used  by  ASAPl.  All  files  are  either  formatted  (F)  or  unformatted  (U) . 


UfliE _ Type _ Description _ _ 

^  ^  Fine  mesh  first  guess  values  for  interpolation.  Contains: 

^^^ijk‘  X  -  (T.U,V)). 

^^ijk-  j-l.J 

where  1-360,  J-181. 

2  Unpacked  OWE  Level  II  data.  See  Norquist  (1984), 
pp.  17-33. 

3  U  Location  information  for  observations.  Contains: 

tN],  dsi^,  Bj^,  points^]  .  n-l,N. 

4  U  Analysis  grid  first  guess  values.  Contains: 

[(Xijk-  k-l,K,  X-(T.U,V)), 

^^ijk-  k-l.Kj,),  i-l,I,  j-l,J 

where  1-61,  J-62. 


Fine  mesh  topography  for  interpolation.  Contains; 
(Z*^j,  i-l.U,  j-l,j 
where  1-360,  J-181. 

Analysis  grid  values.  Contains: 

^^^ijk-  X-(T,U,V)).  i-l,I,  j-l,J 

where  1-61,  J-62. 
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8 

9 

10 


11 


INPUT 

OUTPUT 

Note : 


U  The  sines  of  the  Gaussian  latitudes  (sin  ^ j ] .  j~l.J 

where  J-62. 

U  The  residuals  used  by  the  analysis.  Contains; 

(0BD(n,NU),  NW-1.75).  n-l.N. 

F  Rejected  data.  One  record  per  datum  rejected 

k.  PC-  Pk'  “k- 

<KUk.v^)  I 

U  If  this  is  the  first  analysis  in  a  sequence,  this  file 

contains  the  NMC  EPE  at  the  mandatory  levels  as  a  function 
of  analysis  grid  latitude; 

((E  f-l,L.  X-  (Z.U.V))),  j-l.J,  where  J-62, 

L-12. 

Otherwise,  this  file  contains  the  EAE  for  height  from  the 

previous  analysis:  k-l,K),  i-1,1,  j»l,J, 

aO  ijk 

where  1-61,  J-62. 

U  The  EAE  and  the  analysis  corrections; 

((E  ),  k-l,K,  X  -  {Z,U,V)), 

a  ijk 

^*ijk'  *  "  {z.u,vJ)J,  1-1,1,  j-l,J 

where  1-61,  J-62. 

F  Standard  input  device.  Contains  namelist  DATIHOP. 

F  Standard  output  device.  Contains  printed  output. 


In  this  table  the  contents  of  each  record  is  described  within  square 
brackets  ( (  ] ) . 


Rev  1.10 


A-62 


Table  U  .1 


I/O  activity.  All  I/O  activity  is  listed  here  except  PRINTS  to  unit  OUTPUT. 


Unit 

Activity 

Occurrences 

1 

Rewind 

SETFG  (1) 

Read 

SETFG  (2) 

2 

Open 

ASAPl  (3) 

Read 

ASAPl  (3).  (6).  <7),  (8) 
MASTORl  (2).  (7) 

MASTOR2  (2).  (7),  (22) 
MAST0R4  (2).  (8).  (9).  ( 
MASTOR6  (3) 

Rewind 

ASAPl  (7) 

Open-Close 

ASAPl  (3).  (6).  (7),  (8) 
MASTORl  (18) 

MASTOR2  (23) 

MASTOR4  (18) 

Close 

MASTOR6  (19) 

3 

Write 

ASAPl  (11) 

4 

Read 

ASAP2  (1) 

5 

Rewind 

SETFG  (1) 

Read 

SETFG  (2) 

FG  (2) 

6 

Write 

ASAP2  (49) 

7 

Read 

ASAPl  (13) 

8 

Write 

ASAPl  (11) 

9 

Write 

REJECT  (2) 

10 

Read 

ASAPl  (14) 

ASAP2  (4) 

11 

Write 

ASAP2  (49) 

Input 

Read 

ASAPl  (2) 
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5. 


Methods 


In  this  section  we  document  certain  methods  which  are  used  in  ASAPl.  This 
section  of  documentation  is  independent  of  the  actual  FX)RTRAN  code.  Some  of 
these  methods  are  used  repeatedly  in  the  code. 

5.1  Conversions 

5.1.1  Vector  wind  to  components 

Given  wind  speed  |V|  and  direction  6  the  eastward  and  northward  wind 
components  are  given  by 

u  -  -  |V|  sin  e 

V  -  -  |V|  cos  0 

Note  that  0  is  the  direction  the  wind  comes  from.  (See  also  Norquist  1986b, 

Eq.  23.) 

5.1.2  Heights  to  temperatures  hydrostatically  (and  vice  versa) 

The  Integrated  form  of  the  hydrostatic  relationship  is 

-  (g/R)  AZ  -  /  T  din  p 

Pi 

where  g  is  the  acceleration  of  gravity,  R  is  the  gas  constant  for  dry  air,  AZ  is 
the  thickness  of  the  layer  (Z^-Z^)  between  pressure  levels  and  p^ 

(Pjj  >  p^) ,  T  is  the  temperature  and  p  is  pressure.  Strictly  T  is  virtual 
temperature.  If  T  is  linear  in  In  p,  then 

•  (g/R)  AZ  -  T  Ain  p 

where  T  is  layer  mean  temperature  is  given  by  (T^  +  T^)/2.  A  similar  equation  is 
given  by  Norquist,  (1986b,  Eq.  5). 

5.1.3  Layer  to  level  temperatures 

A  least  squares  procedure  for  converting  between  layer  and  level 
temperatures  due  to  Flattery  is  described  in  detail  by  Norquist  (1986b,  Appendix 
A,  pp  55-58.)  As  an  alternative,  if  one  level  temperature  is  known. 
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then  all  other  level  temperatures  may  be  inferred  from  a  complete  set  of  layer 
temperatures  if  some  assumption  is  made  about  the  functional  form  of  T(p)  between 
the  levels.  If  T  is  linear  in  in  p  then  one  obtains  the  method  described  by 
Norquist  (1986b,  Appendix  D,  pp  64-65).  In  brief,  if  T  -  a  +  b  in  p  between  p^ 
and  Py,  then  since  Jinx  dx  -  xinx  -  x  we  have 

Pi  Pi 

T  /dp  -  /  (a  +  b  in  p)  dp 

P  P 

•^u  u 

or 

(p^-  Py)T  -  p^(a  +  b(in  '  Py(®  +  P^'^)) 

T  -  a  +  b  in  p 

T,-a  +  b  inp. 
i  ^i 

which  are  three  equations  in  the  three  unknowns  a,  b  and  either  or  . 

E.g. ,  if  Ty  is  known,  then 

a  -  T  -  b  in  p 
u  u 

b  -  (T  -  (p^-  Pj/^Pu  Pil'*"<Pi/Pu>-ll> 


5.1.4  Humidity 

Dew  point,  T^,  dew  point  depression  T-T^j,  specific  humidity  q  and 
relative  humidity  r  are  all  used  in  the  analysis.  These  are  related  by  the 
following  expressions.  (Dutton,  1976,  Chap.  8): 

r  -  c/e^ 

,^iv  ,1  1.. 

es  -  *0  ^’'P  ^r~ 

V  o 

-  6,11  expll9.927A-5443.3618/T) 

®  ”  ®s  ^^D^ 

W  -  (l/«)  e/(p-e);  e  -  tWp/d+eW) 
q  -  W/(l+W);  W  -  q/(l-q) 
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where 

f  -  (molecular  weight  air)/(molecular  weight  water)  -  1.61, 

e  is  the  vapor  pressure,  e^  the  saturation  pressure,  and  W  the  mixing  ratio  (i.e. 
Che  ratio  of  Che  mass  of  vapor  to  the  mass  of  dry  air).  Similar  equations  are 
given  by  NorquisC  (1986b,  Eq.  1-4  and  22)  except  that  volume  mixing  ratio  is  used 
and  the  definition  (  is  reversed  in  Eq.  1-4. 

3.2  Generic  quality  marks 

We  dlstinquish  (a)  quality  indicators  Qj ,  which  are  read  from  the  GWE  Level 
II  data,  (b)  generic  quality  marks,  Q,  which  are  used  in  the  analysis  selection 
procedures  and  (c)  quality  levels,  Qj^,  which  are  used  in  the  buddy  check 
procedures . 

The  quality  Indicators  Qj  are  described  by  UMO  (1978),  e.g.  Table  IV  of 
Appendix  A.  Ordinarily,  Qj  -  0  indicates  no  quality  control  check  (QCC)  was  made 
and  Qj  -  1,2...  indicate  decreasing  quality  was  determined  during  the  QCC. 

However  Qj  has  different  meanings  for  different  data.  In  ASAPl  all  Qj  are 
translated  to  Q,  which  always  have  the  same  meaning  (see  Table  5.1).  Note  if 
Q(Zj)  <  Q(Zj)  then  Zj  is  better  Chan  Zy 
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Generic  Quality  Marks 


Quality 

Mark 

Upper  Air 
Quality 
Indicator 

Heanine 

0 

1 

Value  found  correct  during  QCC. 

1 

0 

No  QCC  was  made. 

2 

4 

Value  found  erroneous  during  QCC,  reconstituted 

value  inserted. 

3 

6 

Value  found  missing  during  QCC,  reconstituted  value 

inserted. 

4 

7 

Value  found  missing  during  QCC,  new  value  assigned. 

5 

5 

QCC  made,  new  value  probably  entered. 

6 

2 

Value  found  suspect  during  QCC. 

7 

8 

Value  found  erroneous  during  limits  checks. 

8 

3 

Value  found  erroneous  during  QCC. 

9 

9 

Value  missing. 

The  quality  levels  used  in  the  buddy  check  are  given  by  Norquist  (1986b, 
PP  7-8). 

{ 

5.3  Buddy  Check  Boxes 

The  following  procedure  determines  the  buddy  check  box  index  B^^,  given  the 
latitude  ^  and  longitude  X  of  observation  (Norquist,  1986b,  p.5): 

(1)  Determine  the  latitude  band  i  -  1,...  18  from  SP  -*  NP  in  10*  increments 
by  i  -  (^  +  90)/10  +  1.  The  center  of  the  latitude  band  is  then 

-  lOi  -  95. 

(2)  Determine  the  number  of  boxes  for  this  band  as  n.  -  36  *  cos  ^  . 

A  c 


0+1 

3(0+3)  +  1 
3(0+3)  +  2 
3(0+3)  +  3 


if  type  - 


6a 
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(3)  Determine  the  longtitude  index. 


j  -■  A/AA  +  1  where  AA  -  360/n^. 

(A)  Set  -  lOOi  +  J. 

Note  that,  for  the  CRAY-1,  A  and  ^  must  be  incremented  by  0.001  to  get 
reproducible  results.  This  should  have  no  effect  since  A  and  ^  are  coded  to 
hundredths  only . 

5.4  Vertical  grid  and  interpolation 

A  sigma  coordinate  system  is  used.  The  interface  value  oA  are  given  by 

pA/p  for  k-1,  K+1.  The  layer  thicknesses  Aa  -  oa  -  oa  are  defined  as  0.075, 
k  ^  k  k  k+ 1 

0.125,  0.150,  0.150,  0.125,  0.075,  6*0.050.  In  addition  an  assumed  underground 

layer  has  Ao^  -  0.075.  Indexing  starts  with  the  ground  k-1  and  ends  with  the  top 

of  the  atmosphere  k  -  K+1.  The  layer  values  are  given  by 

K+1  K+1  1/k 

where  k  -  R/Cp  (Brenner  et  al.,  1982,  Eq.  32;  Norquist,  1986b,  Fig.  1). 

Note  that  the  vertical  indexing  for  height  is  generally  offset  by  1. 

Thus,  Z|V  is  the  height  at 

Vertical  interpolation  is  generall  linear  in  inp;  with  data  and  at  p^ 
and  p^j  we  have 

X-X^  in(p/p^) 

X  -X,  "  in(p  /p.) 

In  ASAPl,  this  relationship  is  often  put  in  point- slope  form  using  the  midpoint 
of  the  interval, 

X-(X^+  X^)/2  +  ((X^-  X^)/fn(p^/p^)]  (In  p  -  1/2  In  p^p^) 
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5.5  Polar  stereographic  projection 

For  gridpoints  with  1^  I  S  70*  a  polar  stereographic  projection  is  used  for 

S 

data  selection,  buddy  checking  and  performing  the  analysis.  Map  factors  are  used 
throughout  following  Bergman  (1979)  although  Dey  and  Morone  (1985)  suggrci  that 
for  a  polar  cap  the  map  factors  are  not  necessary. 

The  conversion  of  (A,^)  coordinates  to  (x,y)  coordinates,  in  a  system  in 
which  the  positive  x  direction  points  from  the  pole  towards  Greenwich,  is  given 
by  (Bergman,  1979,  Eqs.  3. 3-3. 4) 

m  -  2/(1  +  sin  (±  4)) ■ 

X  -  am  cos  ^  cos  A 


y  -  ±  am  cos  4  A 


where  m  is  the  map  factor,  a  is  the  earth's  radius  and  the  +  (or  -)  corresponds 
to  the  Northern  (or  Southern)  Hemisphere.  The  distance  d^j  between  two  points 
indexed  i  and  J  is  then  given  by  (Bergman,  1979,  Eq.  3.5) 


where 


2  1  2  2 

**11  “  -  ^^*i*  ’‘i'  *  '^i'  ^ 

J  m  ■'  ■’ 

m  -  (nij^  +  mj)/2. 


In  the  coordinate  system  chosen  for  the  projection,  only  the  winds  at  90*W 
are  unchanged  from  the  (A,^)  system.  The  others  must  be  rotated  by  (Norquist, 
1986b,  Eq.  6) 


u' 

slnA 

±cosA 

'  • 

u 

V' 

TcosA 

sinX 

4 

V 

where  (u*  ,v' )  is  the  wind  in  the  stereographic  projection  and  (u,v)  is  the  wind 
in  the  (A,^)  system.  Again  +  (or  -)  corresponds  to  the  Northern  (or  Southern) 
hemisphere.  The  reverse  conversion  is  given  by  the  transpose  of  the  above 
rotation  matrix: 
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I  sinA  TcosA 
[±cosA  sinA 

5.6  Coefficient  of  geostrophy 

The  coefficient  of  geostrophy  G  decouples  the  wind  and  height  analysis  in 
the  tropics.  Bergman  (1979,  Eq .  3.20)  uses 

G  -  1  -  exp  (1*1/20) 

while  Dey  and  Morone  (1985,  Eq.  BIO)  use 

0  0<|*1<10 

G  -  •  (cos  (12|*(-300)  +  l]/2  10<|*|<25 

^  25<|*|<90 

ASAPl  uses  (Norquist,  1986b,  pl4,  and  Fig.  2) 

(*)  G  (25)  0sl*l<25 

A 

^  25<|*|<90 
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